Direct numerical simulations for hybrid rocket boundary layers: Performance modeling and scaling

Abstract

This paper presents a comprehensive performance and scaling analysis of direct numerical simulations for reacting boundary layers, focusing on slab burner configurations. Using a PETSc-based finite volume CFD framework, the study evaluates the scalability and computational cost of flow, chemistry, and radiation evaluations across 2D and 3D simulations. Polymethyl methacrylate (PMMA) is the fuel with pure O₂ as the oxidizer, modeled using a detailed chemical kinetics mechanism with 113 species and 660 reactions. A ray-tracing-based radiation solver, designed for distributed memory applications, is implemented to model radiation heat transfer. Parallel scalability is analyzed for the coupled flow, chemistry, and radiation heat transfer processes. Weak and strong scaling studies are conducted on up to 15,000 computational ranks, revealing robust performance when flow cells exceed 200 per rank. Chemistry evaluations dominate the computational cost in large 3D simulations, accounting for approximately 40% of the total runtime, while flow processes contribute around 35%, and radiation solver contributions remain below 10% due to reduced evaluation frequencies. GPU accelerated chemistry evaluation, implemented with Zero-RK, demonstrates significant promise, achieving up to a 4x speedup for workloads exceeding 30,000 cells per GPU. However, diminishing returns are observed for smaller workloads due to CPU-GPU communication overhead. This study identifies key challenges, including memory bottlenecks and the effects of domain partitioning on flow scalability, while highlighting the potential of GPU-accelerated chemistry to reduce computational costs. These findings provide realizable run configurations for 2D, 3D, and GPU-accelerated cases, offering insights for optimizing reactive flow solvers.

Keywords

Performance benchmarking scalable algorithms parallel solvers direct numerical simulation CFD scalability hybrid rockets

1. Introduction

Computational Fluid Dynamics (CFD) solvers have become essential tools across various industries and academic disciplines. CFD simulations employ numerical methods to solve and analyze the behavior of fluid flows for problems where analytical solutions do not exist. Insights into fluid behavior can enable optimization of engineering designs as well as development and validation of empirical models. The complexity and associated computational cost of CFD simulations are highly dependent on the type of model being used. Reynolds Averaged Navier Stokes (RANS) simulations provide time-averaged flow information using sub-scale models to model turbulent fluctuations. Large Eddy Simulations (LES) resolve the larger scale turbulent eddies while modeling smaller-scales with various sub-grid models. Direct Numerical Simulations (DNS), the main interest of this study, resolve all spatial scales and therefore require no sub-grid models (Domingo and Vervisch, 2023). The high spatial and temporal resolutions required for DNS lead to substantial computational resources, typically facilitated by High-Performance Computing (HPC) systems.

Current large HPC systems have reached and even surpassed exascale capabilities (Christoudias et al., 2021; Domingo and Vervisch, 2023; Lee et al., 2013; Xie et al., 2024), primarily by increasing the number of computing cores in their clusters. As core-frequency-rates have plateaued in recent years, the primary way to decrease time to solution in CFD simulations has been through increased parallelization (Fischer, 2015). Parallelism is typically achieved by distributing spatial Degrees of Freedom (DOF) across processors for simultaneous computation. Computational speedup achieved through increased parallelism at fixed problem size is a key metric of parallel efficiency and is also referred to as strong scaling performance. In practice, multi-physics CFD codes fall short of ideal speedup, primarily due to communication overhead between domain partitions. Additionally, for CFD codes, memory bandwidth and latency represent significant bottlenecks for intra-node scalability (Chang et al., 2018; Stone et al., 2018). Once the memory bandwidth limit is reached, inter-node parallel efficiency is usually limited by the communication overhead through the network. Communication overhead generally increases with increased partitioning, as the ratio of local computations to communication decreases. Counteracting this tendency requires efficient domain distribution through sophisticated algorithms and techniques to balance computational load (Hassan et al., 2020; Mubarak et al., 2017; Xu et al., 2016).

The difficulty of achieving high efficiency in numerical simulations is compounded by the diverse nature of HPC hardware. Processing architecture and network interconnect characteristics affect performance, leading to significant benefits from tailored algorithms (Mira et al., 2023). Recent trends in computing architectures, particularly the increasing prominence of Graphics Processing Units (GPUs), have significantly enhanced the feasibility of large-scale CFD simulations by reducing time and cost for a broad range of applications (Huismann et al., 2018; Oyarzun et al., 2020; Xu et al., 2016). However, due to the low-bandwidth/high-latency nature of GPU-CPU communication, Single Instruction Multiple Thread (SIMT) methods are necessary to fully utilize the performance of these architectures, instead of the Single Instruction, Multiple Data (SIMD) methods typically implemented on CPUs (Oyarzun et al., 2020; Stone et al., 2018). Regardless, GPU implementations for reacting flows continue to gather interest in the CFD community (Fischer, 2015; Rao et al., 2022; Stone et al., 2018).

Reactive CFD poses unique challenges due to the complexity of chemical reactions and their interactions with advective and diffusive processes (Domingo and Vervisch, 2023; Guan et al., 2023; Mira et al., 2023; Rao et al., 2022; Wang et al., 2024). Detailed chemical kinetic mechanisms, required to sufficiently describe the pyrolysis and oxidation of complex fuels, add significantly to the number of degrees of freedom and computational cost (Barwey and Raman, 2020; Domingo and Vervisch, 2023; Mira et al., 2023). For reacting DNS, the increase in degrees of freedom per cell results in a decrease in the number of cells from $O (1 0^{10})$ to $O (1 0^{7})$ for realizable simulations, limiting these efforts to smaller physical domains due to DNS resolution requirements (Domingo and Vervisch, 2023). Additionally, reacting systems often introduce numerical stiffness, further increasing the computational cost (Ghioldi and Piscaglia, 2024; McNenly et al., 2015; Stone et al., 2018). A common method for handling stiffness is operator-splitting, where the fluid transport and chemical kinetics are solved independently. Operator-splitting also allows one to take advantage of GPU acceleration, by integrating multiple reacting ODE systems simultaneously (Stone et al., 2018). Load imbalance is also a more acute issue for reactive CFD, as the computational effort required for chemistry integration can vary greatly across the simulation domain (Desai et al. 2021).

Performance measurement and optimization of CFD codes are crucial for achieving efficient and scalable simulations. These studies provide detailed information on the overall behavior of the code, and help identify and eliminate bottlenecks in simulations. Garcia-Gasulla et al. conducted performance studies for three different CFD frameworks: OpenFOAM, Alya, and CHORUS, and identified bottlenecks and provided insights for enhancing performance and scalability, such as optimizing the memory access patterns and improving parallel efficiency (Garcia-Gasulla et al., 2020). Another scalability study by Hassan et al. (2020) focused on the scalability of ANSYS Fluent applications, revealing that Intel MPI outperformed other MPI libraries in terms of core speedup and solver efficiency, emphasizing the importance of optimizing communication strategies to minimize latency and maximize bandwidth utilization (Hassan et al., 2020). Furthermore, Karp et al. analyzed the performance of a CFD solver on modular HPC architectures, emphasizing the need to optimize accelerator usage by efficiently combining CPU and GPU resources to accelerate simulations (Karp et al., 2024).

The primary domain of interest in this study is reacting boundary layer theory and small-scale slab burner simulations that capture the same underlying physics as those governing hybrid rocket motors. Slab burner simulations are challenging due to high-speed reactive flows, thermal radiation, multiphase interactions, and complex boundary conditions required to accurately model the boundary layer (Domingo and Vervisch, 2023). To address these challenges, all simulations in this study are conducted using the reactive finite volume flow solver ABLATE (Ablative Boundary Layers at the Exascale) (CHREST, 2024).

The main goal of this study is to determine the performance and scalability of two- and three-dimensional slab burner simulations using the ABLATE framework. Specifically, the study aims to characterize the performance and scalability of the relevant physics models, optimize resource allocation and balance time to solution in the slab burner simulations. First, a single-node scalability study is conducted to assess solver behavior, and required memory bandwidth for the simulations. Strong and weak scalability of ABLATE is evaluated separately for the flow, chemistry, and radiation solvers in 2D and 3D simulations. Scalability is measured on two different computational platforms, Dane and Lassen, in the Lawrence Livermore National Laboratory (LLNL) Open Computing Facility. On Lassen, the GPU accelerators are used for chemistry evaluation through the Zero-RK chemical integrator package (McNenly et al., 2015). The overall speedup achieved through GPU-acceleration is analyzed, considering trade-offs between the complexity of the chemistry, flow, and radiation solvers. The novelty of this study is in the comprehensive performance analysis of a multi-physics CFD code, including high-speed flows, reactive thermochemistry, and thermal radiation, as applied to the slab burner problem. While previous studies have explored the scalability of high-speed flows (AlOnazi et al., 2015; Chang et al., 2018; Räss et al., 2022; Wu et al., 2024; Xie et al., 2024), reactive flows (Guan et al., 2024; Henry De Frahan et al., 2022; Karp et al., 2024; Rao et al., 2022; Stone et al., 2018; Wang et al., 2024), and radiation modeling individually (Li et al., 2023; Mylotte et al., 2025) in different contexts, to the author’s knowledge no one has examined all three simultaneously.

The rest of the study is organized as follows. Section 2 presents the governing equations solved in ABLATE. Section 3 details the numerical techniques implemented in ABLATE to solve flow, finite-rate chemistry, thermal radiation. Section 4 details the performance metrics and the overall scaling strategy used in this work. Section 5 describes the physical problem setup, and gives representative results from slab burner simulations. Section 6 presents and analyzes performance and scaling results for different slab burner configurations. Section 7 summarizes the conclusions and future work.

2. Mathematical formulation

2.1. Flow

The gas phase system for the slab burner environment is described by the compressible Navier-Stokes equations for a reacting flow, coupled with thermal radiation heat transfer,

\begin{aligned} \frac{\partial ρ}{\partial t} + \frac{(\partial ρ u_{i})}{\partial x_{i}} = 0 \\ \frac{\partial ρ u_{j}}{\partial t} + \frac{\partial (ρ u_{i} u_{j})}{\partial x_{i}} = - \frac{\partial p}{\partial x_{j}} + \frac{\partial τ_{i j}}{\partial x_{i}} \\ \frac{\partial (ρ Y_{k})}{\partial t} + \frac{\partial (ρ u_{i} Y_{k})}{\partial x_{i}} = \frac{\partial}{\partial x_{i}} (ρ D_{k} \frac{\partial Y_{k}}{\partial x_{i}}) + {\dot{m}}_{k}^{'''} \\ \frac{\partial (ρ e_{t})}{\partial t} + \frac{\partial (ρ u_{i} H_{t})}{\partial x_{i}} = \frac{\partial}{\partial x_{i}} (u_{j} τ_{i j}) + \\ \frac{\partial}{\partial x_{i}} [(k \frac{\partial T}{\partial x_{i}} + \sum_{k = 1}^{N} ρ D_{k} \frac{\partial Y_{k}}{\partial x_{i}} h_{k})] \\ + \sum_{k = 1}^{N} {\dot{m}}_{k}^{'''} h_{f, k}^{°} - \frac{\partial}{\partial x_{i}} (ρ Y_{C (s)} V_{T, i} h_{C (s)}) - \nabla \cdot {\vec{q}}_{r} \end{aligned}

where, k is the thermal conductivity,

{\dot{m}}_{k}^{'''}

is the mass consumption or production rate of the k^th species from chemical reactions,

e_{t} = e + u_{i}^{2} / 2

is the total sensible energy, H_t = e_t + p/ρ is the total enthalpy,

h_{f, k}^{°}

is the heat of formation of the k^th species, and

τ_{i j} = - 2 / 3 μ \partial u_{k} / \partial x_{k} δ_{i j} + μ (\partial u_{i} / \partial x_{j} + \partial u_{j} / \partial x_{i})

is the viscous stress tensor. V_T,x is the thermophoretic velocity associated with soot thermophoresis (Zimmer, 2016). Equal diffusivities are assumed for each species that vary with temperature, and are estimated assuming unity Lewis number, where the temperature-dependent viscosity is estimated using Sutherland’s law.

2.2. Chemistry

The molar production rate for the ith species, ${\dot{ω}}_{i}$ , is defined as,

{\dot{ω}}_{i} = \sum_{j = 1}^{N} ν_{i j} q_{j}

(1)

where ν_ij is the stoichiometric coefficient for species i, and q_j is the net rate of progress for the j^th reaction, defined as

q_{j} = k_{f, j} \prod_{i = 1}^{N} {(\frac{ρ Y_{i}}{M W_{i}})}^{ν_{i j}^{'}} - k_{r, j} \prod_{i = 1}^{N} {(\frac{ρ Y_{i}}{M W_{i}})}^{ν_{i j}^{''}}

(2)

where k_f,j and k_r,j are the forward and reverse rate constants, respectively. The rate constants are functions of the thermodynamic state, with most elementary reactions described by Arrhenius rate expressions that are primarily functions of temperature. Pressure dependence is included in some reactions through, e.g., three-body reactions and fall-off reactions.

Consistent with the finite-volume formulation of the solver, operator-splitting treats the finite-rate kinetics approximated as constant-volume, constant-density reactors within each cell. The resulting governing equations are:

\begin{aligned} \frac{d Y_{i}}{d t} & = \frac{M W_{i}}{ρ} {\dot{ω}}_{i}, \\ \frac{d T}{d t} & = \frac{R_{u} T \sum_{i} ω_{i} - \sum_{i} {\bar{h}}_{i} {\dot{ω}}_{i}}{\sum_{i} [\frac{ρ Y_{i}}{M W_{i}} ({\bar{C}}_{p, i} - R_{u})]} . \end{aligned}

(3)

where R_u is the universal gas constant, and

{\bar{h}}_{i}

and C_p,i represent the molar enthalpy and the specific heat of the i^th species. The chemical kinetics mechanisms used in this study was developed by Amiri (2024) to describe the pyrolysis and oxidation. The mechanism contains 113 species and 660 reactions, including critical pathways for PMMA oxidation, necessary for accurate modeling of the ignition and extinction behavior of the fuel Amiri (2024). Additionally, a semi-empirical, six-step, two-equation soot model based on acetylene is implemented. This soot model adds two DOF per cell, representing the soot mass fraction and the soot particle number density. The model describes soot nucleation, agglomeration, surface growth, and oxidation by OH, O, and O₂ with reaction rates taken from Ref. Aphale and DesJardin (2021).

2.3. Radiation

Thermal radiation is modeled by the radiative transfer equation (RTE), which describes the change in radiation intensity traveling through a medium. The intensity, I, at a point x and in a direction (θ, ϕ), changes due to absorption and emission by the medium. For a participating medium, the radiative transfer equation is expressed as:

\frac{d I (x, θ, ϕ)}{d s} = κ [I_{b} (T (x)) - I (x, θ, ϕ)]

(4)

Here, κ is the Planck-mean absorption coefficient, I_b is the blackbody intensity at temperature T(x), and s is the path length. The local Planck-mean absorption coefficient, κ, is evaluated as the sum of the gas absorption coefficient, κ_gas, and the soot absorption coefficient, κ_soot, accounting for both gas and soot radiation. An empirical grey gas model is used to determine the gas absorption coefficient (Mylotte et al., 2025; Zimmer, 2016), while a temperature- and pressure-dependent empirical model, based on the Rosseland constant is used to evaluate the soot absorption coefficient (Modest, 2013; Mylotte et al., 2025). Under the grey gas assumption, the radiative heat flux source term for the finite volume model is approximated as the divergence of the radiative heat flux,

{\vec{q}}_{r}

, given as:

\nabla \cdot {\vec{q}}_{r} = κ (4 σ T^{4} - G)

(5)

where σ is the Stefan-Boltzmann constant, and G is the irradiation, computed as;

G (x) = \int_{4 π} I (x, θ, ϕ) d Ω = \int_{0}^{2 π} \int_{0}^{π} I (x, θ, ϕ) c o s (θ) d θ d ϕ

(6)

3. Numerical methods

The coupled system of equations is solved with an explicit Runge-Kutta integrator using the PETSc framework (Balay et al., 2024a, 2024b) for unstructured meshes (Knepley and Karpeev, 2009; Lange et al., 2016). The advective and diffusive source terms along with the boundary conditions and the auxiliary variables, such as temperature are computed in every RK step. The small DNS time-step dictated by the CFL condition on the order of $O (1 0^{- 9})$ to $O (1 0^{- 7}) s$ , allows for looser coupling between chemistry, soot and the flow. As a result, the chemistry source terms are evaluated once at the beginning of each time-step and are held constant throughout the step. Similarly, the soot source terms for solid carbon, number density, and oxidation reactions are linearized in the same manner. Small time-steps also allow for relaxed coupling with radiation processes, as changes in the temperature and soot fields happen on slower timescales, allowing the radiation source terms to be held constant over multiple time-steps. Time-averaged radiative heat flux at the slab burner surface was observed to be insensitive to the frequency of radiative evaluations, remaining consistent with as few as one update every 100 time steps, or around once per microsecond. The ABLATE time-step sequence is summarized in Algorithm 1.

3.1. Flow

For the slab burner simulations, the convective fluxes are discretized using an AUSM + UP flux vector splitting method (Chang and Liou, 2007). Advection Upstream Splitting Methods (AUSM) separate the numerical fluxes into convective and pressure components, which allows for accurate capture of shocks and discontinuities in the flow field. Specifically, the AUSM + UP scheme extends the AUSM family schemes to accurately solve flows at all speed regimes. Pressure and velocity corrections in the method help maintain stability and accuracy, particularly in low Mach number flows. Molecular fluxes are estimated using a first-order difference scheme.

3.2. Chemistry

An open-source chemical kinetic solver, Zero-RK (McNenly et al., 2015; McNenly et al), has been implemented in ABLATE to compute chemical source terms. Zero-RK is designed to efficiently integrate large kinetic mechanisms by utilizing adaptive preconditioners within its implicit ODE solvers. Traditional approaches to chemical kinetics rely on dense matrix operations that become computationally expensive for models over tens of species. Zero-RK addresses this challenge, significantly reducing the computational cost without compromising accuracy McNenly et al. (2015). The Zero-RK plug-in for CFD solvers offers streamlined integration into reactive flow solvers as well as multiple time integration methods suited to specific numerical regimes.

3.2.1. Chemistry solver benchmark

A constant-volume, zero-dimensional ignition problem is used as a benchmark to validate the chemical kinetic integrator implementation in ABLATE. The 113-species kinetic mechanism described in Section 2.2, is used for the benchmark problem. The initial gas mixture consists of MMA vapor and O₂ in stoichiometric proportions, initialized at 1400 K. Figure 1(a) shows the time evolution of temperature and major mass fraction profiles during ignition. A constant large time-step of 1 × 10⁻⁶ s is used to highlight potential differences between the integrators, with identical integration tolerances set for both solver. After ignition, the solvers converge to the same equilibrium state. Figure 1(b) shows the relative error between the baseline Cantera (Goodwin et al., 2017) simulation and Zero-RK implementation in ABLATE. Overall, only minor differences are observed in the temperature and mass fraction profiles, validating the implementation. The largest error, around 1.5%, occurs during the ignition step between Zero-RK and Cantera. In practice, however, significantly smaller differences are expected in slab burner simulations when using smaller time-steps.

3.2.2. GPU accelerated chemistry calculations

Chemical kinetics evaluations are the dominant computational cost for large-scale reactive flow simulations, due to the large number of reaction species necessary to predict burner dynamics (Barwey and Raman, 2020). Zero-RK leverages GPU accelerators, increasingly more common in HPC systems, by simultaneously advancing ODEs for many independent reacting cells (Gao et al., 2018; Oyarzun et al., 2020; Xu et al., 2016). The impact of GPU acceleration on simulation cost was evaluated by solving chemistry on the GPU, while flow and radiation source terms were computed on the CPUs.

3.3. Radiation

Computational models of RTE must accurately capture the interactions between radiation and the participating medium while being computationally efficient. ABLATE implements a ray decomposition method to improve scalability by segmenting rays into sub-domains. This enables distribution of the RTE calculations across processors while limiting inter-processor communication. A detailed performance analysis of the radiation solver in ABLATE is provided by Mylotte et al. (2025).

3.4. Discretization

For multi-core applications, the computational domain is partitioned into subdomains using the PETSc DMPlex method (Balay et al., 2024b). Communication between subdomains is handled via MPI Forum (1994) (MVAPICH 2.3.7), in which a number of ghost cells are used to exchange the required information. Since the flux evaluation time is not dependent on the local state of the cells, this partitioning method generally results in good load balancing. The cost associated with chemistry, however, is highly dependent on the local thermodynamic state. Fuel pyrolysis reactions near the fuel surface and oxidation reactions in the flame introduce numerical stiffness, which slows down integration in the affected cells. To address this, Zero-RK implements an active load balancing algorithm for both CPU and GPU implementations. Work is distributed among the MPI ranks based on the overall integration time taken by each rank during the previous flow time-step. This is especially important in GPU-accelerated runs, where performance differences between ranks can be significant. In the Zero-RK GPU implementation, each device is explicitly bound to a single MPI rank, rather than distributed in a round-robin fashion. GPU-bound ranks can evaluate chemical kinetics significantly faster than CPU-only ranks, so the load balancer assigns more reactors to them. To maximize resource utilization, MPI ranks without GPUs do not remain idle, they concurrently process their assigned reactors using CPU cores.

4. Performance metrics

Understanding the performance metrics of large applications, such as DNS solvers, is significantly more complex and challenging (Chang et al., 2018), compared to micro-benchmarks, where detailed analysis of memory footprints, arithmetic intensity and communication patterns are feasible (Chang et al., 2018; Mylotte et al., 2025). Several combined performance metrics exist for large codes, including those based on floating-point operations, memory usage, communication, memory bandwidth, and latency (Abdelsamie et al., 2016; Chang et al., 2018; Li et al., 2023).

The implemented and practical scalability of ABLATE is evaluated using various benchmark simulations. Implemented scalability refers to the scalability of the analytic implementation of the core flow solver, specifically, the linear algebra routines, the basic EOS evaluations, flux calculations and time-stepping. The goal is to measure the performance of identical operations across different processor counts, independent of solver convergence. This is assessed using laminar, low-Reynolds-number, cold flows on a structured box mesh with a perfect gas EOS. These conditions isolate the solver from the effects of preconditioners in the explicit chemistry solver, the iterative temperature solver within the EOS, and uncertainties introduced by the ray decomposition method in the radiation solver. Practical scalability studies incorporate detailed physics models and evaluates scalability under conditions similar to large scale simulations. These studies are conducted using 2D and 3D simulations on unstructured grids, including models with an ideal gas EOS, reactive chemistry, and thermal radiation. Unlike implemented scalability, which isolates core solver operations, practical scalability accounts for the complex interactions between physics components. This introduces additional constraints such as time-step limitations for flow and chemical stability, memory bandwidth and footprint limitations, and practical considerations tied to solution accuracy and robustness. In this work, practical scalability is assessed by segmenting the simulation into three main physics models: (1) flow, (2) chemistry, and (3) radiation, each with its own dominant computational costs.

The computational rate of the calculations, defined as the total DOF solved per unit of time per core is used as a primary metric of comparison for the sub-models (Chang et al., 2018). The DOF for the flow and chemistry solvers, corresponding to equations (1) and (3) respectively, are defined as follows:

\begin{align} {D O F}_{flow} & = N cells \times (N_{species} + 5) \end{align}

(7a)

\begin{align} {D O F}_{chem} & = N cells \times (N_{species} + 1) \end{align}

(7b)

where N_cells is the number of cells in the domain and N_species is the number of species in the chemical mechanism. The computational rate is then computed as:

Rate = \frac{[Total  DOF  solved]}{[Total  time] \times [# cores]}

(8)

The rate provides a good indication of the computational intensity and efficiency of the solver by quantifying how effectively computational resources are utilized to solve the problem. The second key metric is the parallel efficiency (E_p), which indicates how well the solver scales with increasing core counts. This helps identify the point of diminishing returns, where adding more computational resources yields limited performance gains. Parallel efficiency for strong and weak scaling is calculated using equations (9a) and (9b), respectively.

\begin{align} {Efficiency}_{strong} & = \frac{[Speedup for N nodes]}{[N cores]} \end{align}

(9a)

\begin{align} {Efficiency}_{weak} & = \frac{[Time for 1 core]}{[Time for N cores]} \end{align}

(9b)

4.1. Hardware

The performance of the software is tested on two HPC platforms at LLNL: Dane and Lassen. The specifications of the systems are summarized in Table 1.

Table 1.

Compute node specifications for Dane and Lassen.

	Dane	Lassen
CPU architecture	Intel sapphire rapids	IBM Power9
Cores/Node	112	40
GPU/Node	0	4
CPU memory/Node (GiB)	256	256
GPU architecture	NA	Nvidia V100

Each Dane node is equipped with two Intel Xeon Platinum 8480, 56 core CPUs and 256 GB of DDR5 memory. Figure 2 shows the measured memory bandwidth for 1, 2, 3 and 4 nodes of Dane, varying from 1 core to full occupancy using the STREAM benchmark test (McCalpin et al., 1995). Results show that the average memory bandwidth per core is 275 GB/s per node. Each Lassen node has 40 Power 9 CPU cores with 256 GB of memory and 4 V100 GPUs. The reported peak single CPU memory bandwidth of Lassen is 170 GB/s.

5. Physical problem setup

A slab burner model is selected as the main benchmark problem (Retfalvi et al., 2024), shown in Figure 3(a). The fuel is PMMA, assumed to depolymerize into MMA gas at the surface. The extents of the domain are 25.4 mm × 25.4 mm × 200 mm. The dimensions of the fuel slab are 24.5 mm wide, 9 mm high and 60 mm long with a 45° ramp at the leading edge, starting at 20 mm form the inlet. The initial 2D mesh and boundary conditions are shown in Figure 3(a). The base mesh contains 2892 interior flow cells and 324 boundary cells. A mesh refinement study showed that a 3-level refinement of the mesh shown in Figure 3(a) provides sufficient resolution for grid independent results for the wall heat flux and flame attachment location. After the 3-level refinement, the 2D mesh, used for all the strong scaling tests, consists of 185,088 interior flow cells with a maximum mesh size of 0.1 mm at the wall. The 3D mesh is generated by extruding the 2D mesh in the cross-stream direction using a constant grid spacing. The 3D mesh contains approximately 890,000 flow cells, with an approximate mesh size of 0.4 mm at the wall. The authors acknowledge that the dominant Kolmogorov length scales for the slab burner range from $O (1 0^{- 7}) m$ to $O (1 0^{- 6}) m$ , which are significantly lower than $O (1 0^{- 4}) m$ to $O (1 0^{- 3}) m$ mesh sizes used in this study. Consequently, results from slab burner runs are treated as under-resolved DNS results. However, the aim of the paper is to study the scalability of ABLATE’s different sub-processes and their interactions, rather than resolving the turbulent scales. A fully-developed pipe-flow velocity profile is used as the inlet boundary condition with maximum velocity of 6.28 m/s. Pure O₂ is used as the inlet condition at 300 K, and 1 atm pressure. The fuel boundary is assumed to be at the sublimation temperature of PMMA, 653 K, set as adiabatic, such that the incident heat flux from radiative and conductive heat transfer is assumed to entirely be used to vaporize the fuel. At the outlet, a local one-dimensional inviscid (LODI) boundary condition is applied Poinsot and Lelef (1992). No-slip, isothermal boundary conditions are applied for the walls at 300 K. During the initial transient, the fuel ignites due to an added temperature source and develops into a statistically stationary, stabilized flame. Figure 3(b) shows the temperature field at t = 0.0125 s after ignition. This developed flame state t = 0.0125 s is used as the initial condition for all of the practical scaling studies.

6. Results

Computational costs of the flow and radiation source terms are independent of the time-step, but the cost to compute chemical source terms depends on sub-stepping within the chemical ODE solver, which is highly dependent on the local state. To ensure consistency in scaling results, all benchmarks, with the exception of the large milestone (Section 6.5), are run with a constant time-step of 10⁻⁷ s, while the large 3D simulations use an adaptive time-stepper.

6.1. Implemented scalability

The implemented scalability of the flow solver is evaluated through weak and strong scaling tests using rectangular domains with structured meshes under laminar conditions and a perfect gas equation of state. Figure 4(a) shows the strong scaling results using a 1000 × 350 × 350 structured block mesh. Correspondingly, the number of flow cells per MPI rank decreases from approximately 23,500 to about 8500 between the 4-node and 128-node simulations on Dane. Close-to-linear parallel speedup is achieved, indicating minimal communication overhead. The parallel efficiency decreases to 75% for the largest 128-node (14,336-core) simulations. This drop in the efficiency with increasing ranks is expected due to the domain decomposition and the resulting increased communication between sub-domains. Increased discretization requires more ghost cells, which increases the overall memory as well as the communication. Figure 4(b) shows the computational rate for the weak scaling study of the flow solver from 4 to 128 nodes. Each run has 69,754 cells per MPI rank. The computational rate remains constant for all ranks, indicating excellent implemented scalability, and over 95% parallel efficiency achieved.

6.2. Practical scalability

Practical scalability is assessed through 2D slab burner simulations. These simulations of high-speed reacting flows, with detailed chemistry and thermal radiation, accurately capture the flow physics of slab burner experiments (Retfalvi et al., 2024). Figure 5(a) shows the single node parallel speedup, and Figure 5(b) presents the strong scaling parallel efficiency of the 2D slab burner runs. Figure 5(a) shows sub-optimal scaling between 28 and 112 ranks, as the overall speedup diverges from the ideal line. Consequently, the parallel efficiency drops to 0.6 or 60%, as the memory channels on the node become saturated, highlighting the memory constraints on a node as the primary performance limiter. This is attributed to two main reasons. Firstly, as Figure 2 illustrates the available memory bandwidth plateaus after approximately the 20th rank, and simulations become memory bandwidth limited. Secondly, increased discretization raises the overall memory footprint due to the larger number of ghost cells.

The required memory bandwidth for the 112-rank simulation is determined by distributing the processes across multiple nodes while keeping the total number of ranks fixed at 112, thereby increasing the available memory bandwidth. Figure 6 shows the speedup achieved as a function of increasing memory bandwidth for the 112-rank case. The results indicate that all processes experience speedup as bandwidth increases, up to the use of 4 nodes. Beyond this point, no significant performance gains are observed. As shown in the STREAM test in Figure 2, the available memory bandwidth for 4 nodes is approximately 1200 GB/s, suggesting the effective bandwidth requirement for the 112-rank case. This observation is consistent with earlier findings, as this configuration has 28 MPI ranks per node, which is approximately the threshold beyond which single-node simulations also become memory bandwidth limited.

Figure 7(a) shows the parallel speedup for the multi-node simulations. Since simulations are assumed to be memory bandwidth-limited beyond the first node, the strong scaling results primarily estimate the communication overhead. Initially, near-linear performance is observed for all processes. The overall parallel speedup begins to decline around 1000 ranks, at which point each rank has a relatively low work, around 200 flow cells. Figure 7(b) shows the parallel efficiency achieved by the multi-node simulations. Chemistry exhibits the best parallel efficiency, dropping down to only 90% at 1000 ranks and 50% at 10,000 ranks. Local chemical integration doesn’t have any communication overhead, the decreased parallel efficiency is due to active load balancing. The flow and radiation processes exhibit lower parallel efficiency, reaching approximately 20% at 15,000 ranks. For the flow, the low parallel efficiency is primarily due to increased communication overhead. Additionally, the increased discretization of the flame region contributes to load imbalance in portions of the flux calculations that depend on the chemical state. This imbalance is most pronounced in the face projections, where temperature is recovered from the conserved total energy and species fields using a nonlinear root-finding algorithm. To accelerate this EOS solve, the cell temperature from the previous time step is used as the initial guess. The low efficiency for the radiation process is attributed to increased communication costs in larger scale simulations, consistent with prior findings showing the radiation solver to be both latency- and communication-bound (Mylotte et al., 2025).

The weak scalability of the slab burner runs is assessed through five simulations, with varying discretizations, initialized using nearest node interpolation from the same starting frame, using 1 × 10⁻⁷ second as the constant time-step. Each refinement level involves a cell split in every direction, resulting in a fourfold increase in the number of cells for each refinement level in two dimensions. Figure 8(a) shows the computation rate over the 50 time-steps. As the first node fills up a slight drop is observed for both the flow and the chemistry rates as the simulation become bandwidth limited. For the multi-node cases, ideal scaling is observed for both the chemistry and the flow processes. A significant drop in the rate of the chemistry is observed at 64 nodes (7168 ranks). At a level 4 mesh refinement, the time-step required by the CFL constraint is around 4 × 10⁻⁸ s. Larger time-steps result in non-physical states, affecting the performance of the chemistry integrator, creating additional overhead. The total computational rate, defined as [DOF_Chemistry + DOF_Flow]/[(t_{Chemistry time-step} + t_{Flow time-step}) × N_rank], closely follows the chemistry rate, indicating that the time-step cost is dominated by chemistry. Figure 8(b) shows the parallel efficiency of the weak scaling study. Parallel efficiency remains constant for the chemistry and flow processes for the multi-node cases. In contrast, the radiative transfer process exhibits steadily declining efficiency, dropping to just 5%. To investigate this decline, Figure 8(c) shows the load balance, defined as the ratio of time taken on the fastest rank to the time taken on the slowest rank. Load balance remains relatively constant for both the chemistry at 0.91 and the flow at 0.8 across all ranks. However, the load balance for radiation drops to 0.5 at 64 nodes. This imbalance is attributed to the problem geometry, as the mesh is refined, subdomains near the domain corners contain fewer ray segments than those near the center, leading to the imbalance.

Overall, ABLATE demonstrates robust practical strong and weak scalability for the spectrum of interest. The study shows that the optimal configuration for these 2D simulations on Dane only utilizes 28 ranks, at which point the memory channels are not saturated. However, single node DNS simulations are not viable due to the long time to solution. Multi-node simulations show good parallel scalability for simulations when each MPI rank has over 200 cells. This threshold can serve as a guideline for optimizing further simulations to ensure the efficient use of the available computational resources.

6.3. Three-dimensional simulations

Parallel scalability of the 3D simulations is assessed through a strong scaling study. These simulations capture significantly more physical detail, particularly regarding the turbulent boundary layer over the slab, which is necessary for accurate modeling of the reactive boundary layer in high Reynolds number combustion environments. Due to the high computational cost of the 3D ray tracing approach to radiation, and based on results from the 2D simulations, the frequency of the radiative gain evaluations is decreased to a single evaluation every 50 time-steps. For more consistent comparisons with previous runs, the 3D domain is initialized using the fully developed solutions from the 2D simulations. This is done by assigning each cell in the 3D mesh the state of the nearest cell in the 2D mesh, based on cell center distance. Assuming the 2D mesh lies in the x-y plane, this approach effectively extrudes the 2D initial state into the third dimension.

An additional challenge with the practical scalability of the 3D simulations is the overall memory requirement requiring at least 8 nodes to fit the problem. Figure 9(a) shows the strong scaling speedup for the 3D benchmark. Scalability of chemistry is near perfect throughout all ranks, due to the increased workload in the 3D mesh. The flow and the radiation speedup is lower compared to the 2D simulations due to increased overhead associated with additional ghost cells and significantly more rays required for 3D meshes per cell for radiation. Figure 9(b) shows the computational rate for the strong scaling test. As expected, the computational rate decreases for all cases due to the overhead associated with active load balancing in the chemistry and due to the presence of ghost cells for the flow solver. The decline in the 3D chemistry rate occurs later than in the 2D simulations due to the increased workload per MPI rank. Consistent with the rate definition, 3D simulations exhibit lower rates due to the increased number of faces per cell.

6.4. GPU accelerated two-dimensional simulations

The performance of the 2D simulation with GPU-accelerated chemistry evaluations is examined on Lassen, focusing on comparisons with the representative CPU-only cases. As with the CPU cases, data is collected over 50 time-steps, using the same time-steps. All Lassen tests are carried out using the same 2D mesh with the same initial conditions as the 2D CPU runs, representative of the fully developed flame solutions.

Figure 10(a) shows the strong scaling speedup on Lassen. The radiation and flow solvers demonstrate similar speedup trends compared to CPU-only cases. The scalability of chemistry can be divided into three distinct regions. In the first region on the first node, as the number of ranks increases, the number of GPUs assigned remains constant, consequently no significant speedup is expected. Minimal speedup is achieved due to the active load balancing paradigm, as part of the work is assigned to CPUs not bound to GPUs. Linear scaling is observed in the second region for chemistry, as the number of GPUs increases linearly with the number of ranks. In the third region, as the workload per GPU decreases, the communication overhead between GPUs and CPUs begins to degrade the scalability. Figure 10(b) shows the computational rate for strong scaling. The computational rate as expected for the strong scaling test decreases with increasing ranks. For lower ranks, the chemistry rate is significantly higher in the GPU cases compared to the CPU cases, attributed to the relatively large workload given to the GPU. As the number of ranks increases and the problem is further partitioned in a strong scaling sense, the chemistry rate drops sharply. The decline is due to the reduced workload per GPU and the increased communication overhead between CPUs and GPUs. The GPU chemistry rate falls below the CPU rate at approximately 30,000 cells per GPU in these simulations. In contrast, the flow evaluation rate is consistently lower in the mixed CPU-GPU cases compared to the CPU-only runs. This can be explained by the relatively lower bandwidth and computational rate of the CPUs on the Lassen system. This outcome aligns with the design philosophy of HPC-GPU systems like Lassen, where peak computational performance is prioritized for the GPUs, often resulting in reduced memory bandwidth and computational capacity on the CPUs.

Figure 11 shows the speedup per time-step, defined as the ratio of time needed for single time-step on Lassen (GPU) to the time taken for a single time-step on Dane (CPU only). The radiative gain is updated only every 50 flow time-steps. This overall speedup and specifically the chemistry speedup decreases with increasing ranks, while the flow and radiation remains constant at 0.5. Similar to the results shown in Figure 10(b), using 4 ranks where each rank is bound to a GPU, a significant 4-times speedup is observed in the chemistry evaluations, and an overall 1.9-times speedup is observed in the overall time-step. This translates to overall faster time-to-solutions for the smaller simulations. However, as the chemistry rate decreases with increasing ranks, the overall time-step speedup falls below 1 after approximately 20 ranks, indicating that beyond this point, Lassen is slower than Dane.

GPUs have been shown to enhance the performance of ABLATE, but fully leveraging their capabilities requires careful consideration of the problem setup. The goal is to maximize the chemistry workload to offset the lower performance of CPUs running radiation and flow. In larger simulations, this is achieved by increasing the number of cells per rank. Additionally, complex fuels such as paraffin wax, often used in ablating hybrid rockets, lack simple chemical mechanisms Amiri (2024). The chemistry workload for simulations involving these fuels is significantly higher compared to PMMA, making them more attractive for GPU-accelerated simulations.

6.5. Large 3D slab burner simulations

To accurately study the boundary layer dynamics and the fuel response, and to compare to experiments, simulation data is time averaged over multiple flow-through times. Generally, the slab burner simulations advance 0.8 to 1 flow-through per day. In this section, a 3D slab burner simulation, configured to represent a milestone simulation, is evaluated over a single 24-h allocation period. This setup uses the same 3D mesh as previously introduced. The simulation on 16 nodes (1792 ranks) has approximately 500 cells per rank, exceeding the 200 cells per rank threshold, is representative of large milestone setups where simulations are expected to scale and perform most efficiently. Simulation progressed 2.37 ms, equivalent to approximately 1.2 flow-through times, while completing 55,200 time-steps, averaging around 4.2 × 10⁻⁸ s per time-step. Figure 12 shows the relative contributions of each process. The simulation outputs data every 1 × 10⁻⁴ s, and the radiative gains are only evaluated every 50 time-steps. The highest contribution comes from the chemistry process, accounting for approximately 40% of the overall time, while the flow process accounts for 25%, and radiation contributes only 5% due to the less frequent flow evaluations. Results are encouraging, especially considering the complexity of the chemical kinetic mechanism, which includes over 100 species and more than 650 reactions. Solving for the temperature field as part of the auxiliary vector, accounting for 11% of the total allocation time is also a factor in the cost of the simulation. Boundary conditions account for 9% of the overall time. Parallel I/O operations, amount to around 1% of the total simulation time for 15 outputs. Overall, this is considered very acceptable I/O performance, considering each output file exceeds 4 GB.

7. Conclusions

This study provides an in-depth performance and scalability analysis of the ABLATE DNS framework for reacting boundary layers for eventual application to hybrid rockets. The unstructured, high-speed flow solver is fully coupled with the chemistry solver and a detailed ray-tracing-based thermal radiation solver. Chemistry evaluations are conducted using Zero-RK, using detailed chemical kinetics mechanisms with 113 species and more than 660 reactions. The primary goals of this study are to assess the performance of the flow solver, chemistry solver, and radiation solver, and to determine the ideal configuration for efficient resource allocation. The secondary objective is to highlight additional simulation constraints associated with the physics of the problem, providing insights into the challenges of simulating hybrid rocket DNS.

Overall, ABLATE demonstrates robust practical strong and weak scalability for the flow, chemistry and radiation processes. Single-node results show that simulations quickly become memory-limited, as expected for large solvers. The required memory bandwidth is approximately four times greater than the available bandwidth on a single node, resulting in a memory constraint. The communication overhead, characterized through a strong scaling study, shows that simulations scale well up to around 200 cells per rank. The most ideal running configuration for these 2D simulations on Dane only utilizes 28 ranks, at which point the memory channels are not saturated. However, single-node DNS simulations are not viable due to the high time to solution, thus running configurations up to 200 cells per rank are recommended.

GPU-accelerated chemistry evaluations show significant promise, with speedups of up to 4× for smaller simulations. To optimize running configurations for the GPU simulation, the goal is to maximize the chemistry workload to offset the lower performance of CPUs running radiation and flow. For larger simulations, this is achieved by increasing the number of cells per rank. For larger simulations in the strong scaling study, the benefits of GPU acceleration diminish due to the communication overhead between CPUs and GPUs. It is important to note, however, that the chemistry process continues to benefit from GPU acceleration in the earlier stages of scaling.

Future work will focus on further improving the performance of the flow solver and exploring different problem spaces beyond the slab burner, such as DNS of bench-scale hybrid rocket motors and LES of full-scale hybrid rocket engines.

Footnotes

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the United States Department of Energy’s (DoE) National Nuclear Security Administration (NNSA) under the Predictive Science Academic Alliance Program III (PSAAP III) at the University at Buffalo, under contract number DE-NA0003961. The work at LLNL was performed under the auspices of the U.S. Department of Energy by LLNL under Contract DE-AC52-07NA27344.

ORCID iD

Kolos Retfalvi

Author biographies

Kolos Retfalvi is a final-year Ph.D. student in the Department of Mechanical and Aerospace Engineering at the University at Buffalo. His research focuses on high-performance scalability and performance modeling of CFD solvers, chemical kinetics, and uncertainty quantification. He is a member of the Energy Transport Laboratory and Center for Hybrid Rocket Exascale Simulation Technology at the University at Buffalo.

Russell Whitesides has pursued a variety of topics in mechanical engineering R&D in the course of his academic and research career. His doctoral thesis focused on the atomistic chemical mechanisms of soot particle growth in combustion environments. Since joining LLNL he has applied his theoretical and applied knowledge of chemical kinetics and scientific computing platforms towards internal combustion engine simulations with the goal of highly-efficient, clean-combustion for transportation. Dr. Whitesides has also worked with many other groups at LLNL to advance high-performance computing for scientific and engineering simulations.

Matthew G. Knepley is a Professor in the Computer Science and Engineering Department of the University at Buffalo. He received his Ph.D. in Computer Science from Purdue University in 2000. His research focuses on scientific computation, including scalable algorithms and parallel computing, numerical analysis, software development, and computational geophysics, plasma physics, and biology. He is an author of the widely used PETSc library for scientific computing. He was a J.˜T. Oden Faculty Research Fellow at the Institute for Computation Engineering and Sciences, UT Austin, in 2008, won the R&D 100 Award in 2009, and the SIAM/ACM Prize in Computational Science and Engineering in 2015 as part of the PETSc team.

Paul E DesJardin is a Professor at the State University of New York at Buffalo (SUNY Buffalo) with a research focus on experimental and computational studies of turbulent reacting interfaces with applications to propulsion, fire science, biomass energy utilization and hydrogen combustion. He currently serves as Director of the Combustion and Energy Transport Laboratory (CET - https://www.cet-lab.org) and a DOE PSAAP Center for Hybrid Rocket Exascale Simulation Technology (CHREST - https://www.buffalo.edu/chrest.html). He is a Fellow of ASME and an Associate Fellow of AIAA.

Matthew T McGurn has focused his career upon the development of complex simulation frameworks and required interfaces to allow these powerful tools to be used by a wide variety of people. Holding his PhD in Mechanical Engineering (2013) from the University at Buffalo, he is well versed in the mathematics, physics, and software programming methods associated with CFD and FE based formulations.

References

Abdelsamie

Fru

Oster

, et al. (2016) Towards direct numerical simulations of low-mach number turbulent reacting and two-phase flows using immersed boundaries. Computers & Fluids 131: 123–141.

AlOnazi

Keyes

Lastovetsky

, et al. (2015) Design and optimization of openfoam-based cfd applications for hybrid and heterogeneous hpc platforms. ArXiv Preprint arXiv:1505.07630.

Amiri

(2024) Automated Generation of Chemical Kinetic Reaction Mechanisms for Combustion of Hybrid Rocket Fuels - Large Alkanes and Poly Methyl Methacrylate (PMMA). PhD Thesis, The State University of New York, Buffalo.

Aphale

DesJardin

(2021) Two-color pyrometry based flame to fuel surface radiative heat flux diagnostic using flamelets. Combustion and Flame 229: 111395.

Balay

Abhyankar

Adams

, et al. (2024a) PETSc/TAO users manual. Argonne National Laboratory, Technical Report ANL-21/39 - Revision 3.22. DOI: 10.2172/2205494.

Balay

Abhyankar

Adams

, et al. (2024b) PETSc web page. https://petsc.org/

Barwey

Raman

(2020) A neural network inspired formulation of chemical kinetics. ArXiv Preprint arXiv:2008.08483.

Chang

Liou

(2007) A robust and accurate approach to computing the compressible multiphase flow: stratified flow model and ausm+-up scheme. Journal of Computational Physics 225: 840–873.

Chang

Fabien

Knepley

, et al. (2018) Comparative study of finite element methods using the time-accuracy-size (tas) spectrum analysis. ArXiv Preprint arXiv:1802.07832.

10.

CHREST UaB (2024) Ablative boundary layers at the exascale. https://github.com/UBCHREST/ablate

11.

Christoudias

Kirfel

Kerkweg

, et al. (2021) Gpu optimizations for atmospheric chemical kinetics. In: The International Conference on High Performance Computing in Asia-Pacific Region, HPC Asia 2021, Jeju, South Korea, 20–22 January 2021, pp. 136–138. ACM.

12.

Desai

Kim

Song

, et al. (2021) Direct numerical simulations of turbulent reacting flows with shock waves and stiff chemistry using many-core/gpu acceleration. Computers & Fluids 215: 104787.

13.

Domingo

Vervisch

(2023) Recent developments in dns of turbulent combustion. Proceedings of the Combustion Institute 39(2): 2055–2076.

14.

Fischer

(2015) Scaling limits for pde-based simulation (invited). In: 22nd AIAA Computational Fluid Dynamics Conference, Dallas, TX, 22–26 June 2015. American Institute of Aeronautics and Astronautics.

15.

Gao

Grover

Gopalakrishnan

, et al. (2018) Steady-state calibration of a diesel engine in computational fluid dynamics using a graphical processing unit-based chemistry solver. Journal of Engineering for Gas Turbines & Power 140(10): 102802.

16.

Garcia-Gasulla

Banchelli

Peiro

, et al. (2020) A generic performance analysis technique applied to different cfd methods for hpc. International Journal of Computational Fluid Dynamics 34(7–8): 508–528.

17.

Ghioldi

Piscaglia

(2024) A hybrid cpu-gpu paradigm to accelerate reactive cfd simulations. International Journal for Numerical Methods in Fluids 96(8): 1461–1488.

18.

Goodwin

Moffat

Speth

(2017) Cantera: an object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes. Version 2.3.0. DOI: 10.5281/zenodo.170284.

19.

Guan

Abdelsamie

Chi

, et al. (2023) A dedicated reduced kinetic model for ammonia/dimethyl-ether turbulent premixed flames. Combustion and Flame 257: 113002.

20.

Guan

Chi

Liang

, et al. (2024) Revisiting performance of reactivity stratification with hydrogen addition for ammonia combustion. Proceedings of the Combustion Institute 40(1–4): 105381.

21.

Hassan

Rabhy

Mohamed

(2020) Performance study of cfd pressure-based solver on hpc. In: 2020 IEEE International Symposium on Signal Processing and Information Technology (ISSPIT), Louisville, KY, USA, 09–11 December 2020. IEEE.

22.

Henry de Frahan

Rood

Day

, et al. (2022) PeleC: an adaptive mesh refinement solver for compressible reacting flows. The International Journal of High Performance Computing Applications 37(2): 115–131.

23.

Huismann

Lieber

Stiller

, et al. (2018) Load Balancing for CPU-GPU Coupling in Computational Fluid Dynamics. Springer International Publishing, pp. 337–347.

24.

Karp

Suarez

Meinke

, et al. (2024) Experience and analysis of scalable high-fidelity computational fluid dynamics on modular supercomputing architectures. The International Journal of High Performance Computing Applications 39: 329–344.

25.

Knepley

Karpeev

(2009) Mesh algorithms for PDE with Sieve I: mesh distribution. Scientific Programming 17(3): 215–230.

26.

Lange

Mitchell

Knepley

, et al. (2016) Efficient mesh management in Firedrake using PETSc-DMPlex. SIAM Journal on Scientific Computing 38(5): S143–S155.

27.

Lee

Malaya

Moser

(2013) Petascale direct numerical simulation of turbulent channel flow on up to 786k cores. In: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, SC13, Denver, CO, USA, 17–22 November 2013, pp. 1–11. ACM.

28.

Liu

Zhang

, et al. (2023) An efficient heterogeneous parallel algorithm of the 3d moc for multizone heterogeneous systems. Computer Physics Communications 292: 108806.

29.

McCalpin

Whitesides

Flowers

, et al. (1995) Memory bandwidth and machine balance in current high performance computers. IEEE computer society technical committee on computer architecture (TCCA) newsletter. Vol. 2, No. (19-25).

30.

McNenly

Whitesides

Lapointe

Zero-rk: zero order reaction kinetics. https://github.com/LLNL/zero-rk/

31.

McNenly

Whitesides

Flowers

(2015) Faster solvers for large kinetic mechanisms using adaptive preconditioners. Proceedings of the Combustion Institute 35(1): 581–587.

32.

Mira

Pérez-Sánchez

Borrell

, et al. (2023) Hpc-enabling technologies for high-fidelity combustion simulations. Proceedings of the Combustion Institute 39(4): 5091–5125.

33.

Modest

(2013) Radiative Heat Transfer. McGraw-Hill.

34.

MPI Forum (1994) Mpi: A message-passing Interface Standard. Technical report. USA.

35.

Mubarak

Carothers

Ross

, et al. (2017) Enabling parallel simulation of large-scale hpc network systems. IEEE Transactions on Parallel and Distributed Systems 28(1): 87–100.

36.

Mylotte

McGurn

Budzinski

, et al. (2025) Ray decomposition radiation transport for high performance computing. Journal of Computational Physics 521: 113567.

37.

Oyarzun

Chalmoukis

Leftheriotis

, et al. (2020) A gpu-based algorithm for efficient les of high Reynolds number flows in heterogeneous cpu/gpu supercomputers. Applied Mathematical Modelling 85: 141–156.

38.

Poinsot

Lelef

(1992) Boundary conditions for direct simulations of compressible viscous flows. Journal of Computational Physics 101(1): 104–129.

39.

Rao

Chen

, et al. (2022) An investigation of hybrid cpu-gpu solvers for supersonic reacting flow simulation with detailed chemical kinetics. Aerospace Science and Technology 126: 107597.

40.

Räss

Utkin

Duretz

, et al. (2022) Assessing the robustness and scalability of the accelerated pseudo-transient method. Geoscientific Model Development 15(14): 5757–5786.

41.

Retfalvi

Ismael

Budzinski

, et al. (2024) Direct numerical simulation and measurements of a solid fuel slab burner experiment. In: AIAA SCITECH 2024 Forum, Orlando, FL, 8–12 January 2024. American Institute of Aeronautics and Astronautics.

42.

Stone

Alferman

Niemeyer

(2018) Accelerating finite-rate chemical kinetics with coprocessors: comparing vectorization methods on gpus, mics, and cpus. Computer Physics Communications 226: 18–29.

43.

Wang

Chi

, et al. (2024) Direct numerical simulations of turbulent premixed cool flames: global and local flame dynamics analysis. Combustion and Flame 270: 113759.

44.

Karzhaubayev

Shen

, et al. (2024) Direct numerical simulation of three-dimensional particle-laden thermal convection using the lattice Boltzmann method. Computers & Fluids 276: 106268.

45.

Xie

Feng

Huang

, et al. (2024) Extreme-scale direct numerical simulation of incompressible turbulence on the heterogeneous many-core system. In: Proceedings of the 29th ACM SIGPLAN Annual Symposium on Principles and Practice of Parallel Programming, PPoPP ’24, London, UK, 25–28 February 2024, pp. 120–132. ACM.

46.

Hong

, et al. (2016) Optimization of high performance computing cluster based on intel mic. In: 2016 2nd IEEE International Conference on Computer and Communications (ICCC), Chengdu, China, 14–17 October 2016. IEEE.

47.

Zimmer

(2016) Numerical study of soot formation in laminar ethylene diffusion flames. PhD Thesis, Universidade Federal do Rio Grande do Sul.