Stabilizing the Maximal Entropy Moment Method for Rarefied Gas Dynamics at Single-Precision

TL;DR

This work stabilizes the maximal entropy moment method (MEM) for rarefied gas dynamics by introducing gauge transformations, a canonical exponential-family formulation, and a modified Newton solver, enabling reliable single-precision GPU simulations of strong normal shocks (up to Mach $M=10$) with $35$ moments. By transforming to an optimal coordinate system, applying log-sum-exp partitioning, and incorporating a Hermite gauge, the method mitigates ill-conditioning and numerical overflow that previously limited MEM. The study also identifies fundamental resolution limits tied to the mean free path and particle realizability, showing MEM can accurately capture near-continuum shock structures with coarser grids than the mean free path. Overall, the paper presents a practical pathway to deploy MEM on modern GPUs for challenging non-equilibrium flows, with broad implications for multiscale gas-dynamics simulations.

Abstract

The maximal entropy moment method (MEM) is systematic solution of the challenging problem: generating extended hydrodynamic equations valid for both dense and rarefied gases. However, simulating MEM suffers from a computational expensive and ill-conditioned maximal entropy problem. It causes numerical overflow and breakdown when the numerical precision is insufficient, especially for flows like high-speed shock waves. It also prevents modern GPUs from accelerating MEM with their enormous single floating-point precision computation power. This paper aims to stabilize MEM, making it possible to simulating very strong normal shock waves on modern GPUs at single precision. We improve the condition number of the maximal entropy problem by proposing gauge transformations, which moves not only flow fields but also hydrodynamic equations into a more optimal coordinate system. We addressed numerical overflow and breakdown in the maximal entropy problem by employing the canonical form of distribution and a modified Newton optimization method. Moreover, we discovered a counter-intuitive phenomenon that over-refined spatial mesh beyond mean free path degrades the stability of MEM. With these techniques, we accomplished single-precision GPU simulations of high speed shock wave up to Mach 10 utilizing 35 moments MEM, while previous methods only achieved Mach 4 on double-precision.

Figures (8)

• Figure 1: Scatter plots show randomly chosen parameters close to an equilibrium state in the canonical and non-canonical forms, with red and blue dots indicating whether numerical overflow happens when computing moments from these parameters at single precision. (a) The scatter plot of $\beta_5$ and $\beta_6$, components of the randomly chosen parameters in the canonical form defined in \ref{['Maximal Likelihood equation for Exponential Family Distributions: exp family']}. None of them encounter numerical overflow in computing moments. (b) The scatter plot of $\alpha_5$ and $\alpha_6$, components of the randomly chosen parameters in the non-canonical form defined in \ref{['Maximal Likelihood equation for Exponential Family Distributions: exp family ori']}. Many encounter numerical overflow except for a small portion in the third quadrant. These two plots show that the canonical form performs better than the non-canonical form in preventing numerical overflow in the neighborhood of the equilibrium.
• Figure 2: This plot compares the performance of Newton's method and the modified Newton's method in optimizing the parameters corresponding to given moments at single precision. The parameter $\beta_4$ is plotted against the corresponding moment $M_4$ while keeping all other moments as constants. The blue line represents the parameter $\beta_4$ optimized using the modified Newton's method, which is in perfect agreement with the reference values computed at double precision shown as dots on the plot. On the other hand, the red line, representing the parameter $\beta_4$ optimized using Newton's method, fails to match the reference values beyond $M_4 = 4.5$. Overall, this plot demonstrated that the modified Newton's method is more accurate and robust in solving parameters from moments compared to the Newton's method.
• Figure 3: These figures show how gauge transformation speeds up optimization at single precision. It compares the iterations needed to optimize parameters with or without gauge transformation. The goal is to find the parameters for the upper stream flow of normal shock at different Mach numbers. Figures (a) and (c) show the iterations for different Mach numbers and tolerance levels (tol=1e-6 and tol=1e-9). Without gauge transformation, optimization is slow and unstable. With Hermite gauge, optimization is fast and stable. Figures (b) and (d) show the log-mean-square error of the parameters at single precision compared to double precision. The error with Hermite gauge is much lower than without gauge transformation. These figures prove that gauge transformation makes optimization more accurate and faster.
• Figure 4: The plots show the normal shock wave structure at Mach 1.2 computed by the MEM and DSMC methods. The MEM matches the reference DSMC (775 cells) in density, velocity, and heat flux profiles (Figures a-c), while using a coarser resolution of 4 times the mean free path. The DSMC (94 cells) deviates from the reference when the cell width is similar to the mean free path (Figure d). The MEM can capture the flow accurately with resolutions larger than the mean free path hence reduces computation costs, which is hard for mesoscopic methods like the DSMC.
• Figure 5: These plots provide the first prediction of Mach 10 shock wave by MEM utilizing 35 moments of the normal shock wave at single precision, which was too ill-conditioned to be solved previously even with double precision. The density profile in (a) matches DSMC, while the temperature profile in (b) has a higher peak and amplitude near the shock. The stress and heat flux profiles in (c) and (d) also agree with DSMC at and behind the shock, but show higher amplitude in front of the shock. In addition, tiny sub-shocks appear at the upstream end of the stress and heat flux profiles as known artifacts of moment equations. These single precision results indicate that our methods greatly enhanced the numerical stability of the MEM equations, which is necessary to simulate flow at such a high Mach number.