Source-linked AI summary
An implicit multiscale solver for steady kinetic equations based on wave-particle decomposition
Chang Liu, Kun Xu
TL;DR
Steady kinetic computation must capture non-equilibrium transport on meshes larger than the mean free path while addressing slowly decaying macroscopic errors. This paper develops a closed wave–particle formulation and residual-based implicit solver, and reports kinetic fidelity with substantially fewer iterations and less wall time near the continuum limit than the conventional iterative scheme.
Problem
Steady kinetic computations must capture non-equilibrium distributions in rarefied regions while recovering fluid transport on meshes larger than the mean free path, where conventional iterations also leave slowly decaying macroscopic errors.
Source: §1 · p. 2Method
The closed WPD system makes the wave depend on conservative variables and couples macroscopic prediction to microscopic correction through a residual-based implicit iteration.
Source: §6 · p. 41Results
Near the continuum limit, IWPD substantially reduces iteration counts and wall time relative to the conventional iterative scheme, while BGK and Shakhov flow results support the analyzed mechanisms.
Source: §6 · p. 41Takeaways & Limitations
The reported flow tests and parameter studies provide a basis for selecting iteration controls and improving inner solvers.
Source: §6 · p. 41The formal fluid-limit analysis is for a fixed mesh with positive horizons; keeping η finite under mesh refinement means the transport matrix cannot be discarded as an O(ε) perturbation in operator norm.
Source: §4.2 · p. 24
Abstract
from arXiv · showhide
Steady kinetic computations require accurate transport on meshes that do not resolve the mean free path and efficient removal of slowly decaying macroscopic errors. We develop a closed wave-particle decomposition (WPD) and a residual-based implicit solver for this problem. The wave accumulates only the Maxwellian along a local characteristic horizon and is determined by the conservative variables; the particle carries the kinetic complement and collision remainder. Under moment compatibility, the closed system is equivalent to the kinetic equation. For the Bhatnagar-Gross-Krook and Shakhov models, a hydrodynamic expansion gives the wave and particle contributions to the Euler and Navier-Stokes fluxes. The implicit iteration couples macroscopic prediction to microscopic correction through a discrete endpoint term and synchronizes the conservative moments after each accepted correction. Its stationary states satisfy the target discrete kinetic equation independently of the iteration controls. Formal analysis gives the smooth fixed-mesh fluid limit. Linear model analysis relates convergence to the horizon-to-relaxation ratio, incomplete macroscopic solves and explicit collision remainders. One- and two-dimensional computations of a normal shock, internal flows and hypersonic external flows test these properties across rarefied and continuum regimes. The results show kinetic fidelity and substantial reductions in iteration counts and wall time relative to the conventional iterative scheme near the continuum limit. Collision-model comparisons and parameter studies support the predicted convergence mechanisms.
1 Introduction
The paper addresses the need to capture rarefied kinetic behavior on under-resolved meshes while reducing slowly decaying macroscopic errors. It introduces a closed wave–particle decomposition and a residual-based implicit iteration, then tests their behavior across flow regimes.
Contributions
- The wave accumulates only the Maxwellian, while the particle retains the kinetic complement and collision remainder.For BGK, this accumulation choice coincides with the alternative described by Liu and Xu [32].Source: §1 · p. 2
- The closed WPD system couples conservative variables to a particle distribution, and its wave and particle contributions recover Euler and Navier–Stokes fluxes.The ratio η = ℓ/τ controls the decomposition.Source: §1 · pp. 2–3
- The residual-based implicit iteration transfers macroscopic prediction through a discrete endpoint term and synchronizes moments after each accepted microscopic correction.Its stationary states satisfy the target discrete kinetic equation independently of iteration controls.Source: §1 · p. 3
- The analysis examines stationary-state consistency, the formal fluid limit, and how horizon scale, incomplete macroscopic solves, and explicit collision remainders relate to convergence.Source: §1 · p. 3
Numerical tests
- Shock, internal-flow, cylinder, and X38 tests examine kinetic fidelity, hydrodynamic behavior, convergence mechanisms, and comparisons with other methods.Cavity results on grids from a published GSIS study provide an additional method comparison.Source: §1 · p. 3
2 Closed wave–particle decomposition
The closed WPD formulation represents the kinetic distribution using a wave determined by conservative variables and a particle component. Under moment compatibility, the resulting system is equivalent to the kinetic equation and apportions macroscopic fluxes across the two components.
2.2 Closed equations and kinetic equivalence
- Under the moment constraint, the component moments sum to the total conservation law, and the closed WPD system is equivalent to the kinetic equation.The figure depicts the horizon split and the changing transport shares near equilibrium.Source: §2 · p. 3, §2.2 · p. 7, Figure 1 · p. 7
2.1 Kinetic equation and wave–particle components
- For BGK and Shakhov models, the hydrodynamic expansion assigns parts of the Euler and relaxation fluxes to the wave and particle components.The relaxation flux carries the viscous stress and a fraction Pr of the heat flux; the correction flux carries the remainder, while BGK has no correction flux.Source: §2.1 · p. 5
- The wave accumulates the Maxwellian along a characteristic horizon, while the particle retains the attenuated incoming distribution and accumulated collision remainder.The decomposition parameter is η = ℓ/τ; the collision model enters through the remainder.Source: §2.1 · pp. 4, 6
2.2 Closed equations and kinetic equivalence
- The closed WPD system evolves conservative variables and a particle equation, with closure supplied by the wave functional of the conservative state.The moment constraint is essential to the steady solution and is enforced by the implicit algorithm's acceptance step.Source: §2 · p. 3, §2.2 · pp. 6–7
Regularity conditions
The section establishes kinetic equivalence under regularity and moment-compatibility conditions, then describes how the WPD components recover hydrodynamic transport. Its expansion applies under stated smoothness and frozen-relaxation-time assumptions, with a correction when relaxation time varies along characteristics.
Regularity conditions
- Under the stated regularity conditions and moment constraint, solutions of the closed WPD system reconstruct kinetic solutions, and vice versa, for admissible horizons.For a unique kinetic solution, compatible data give the same conservative state and total distribution for different horizons.Source: Regularity conditions · pp. 8–9
Theorem 2.2 (Hydrodynamic expansion of the wave and particle
- With frozen relaxation time and locally constant η, the component expansion holds through first order with O(ε^2) remainders, and the total flux is Navier–Stokes for every η.The expansion assumes the stated smoothness conditions; if τ varies along a characteristic, an additional wave term is compensated by the particle.Source: Theorem 2.2 (Hydrodynamic expansion of… · pp. 11–12
2.3 Wave–particle decomposition in the hydrodynamic limit
- For BGK and Shakhov, the particle share β = (1 + η)e^−η of relaxation transport decreases from one to zero and equals one-half at η ≈ 1.68.As η tends to zero the particle carries all three fluxes; as η grows, the wave carries the Euler and relaxation fluxes while the particle retains the correction share.Source: Theorem 2.2 (Hydrodynamic expansion of… · p. 11
Theorem 2.2 (Hydrodynamic expansion of the wave and particle
- At finite η, the particle contains an equilibrium share that must follow the macroscopic prediction, motivating the endpoint correction in the implicit coupling.Source: Theorem 2.2 (Hydrodynamic expansion of… · p. 12
3 Implicit wave–particle iteration
The discrete formulation couples a Maxwellian wave and kinetic particle through compatible fluxes and residuals. IWPD alternates macroscopic prediction and microscopic correction, using an endpoint term and moment synchronization to preserve the discrete kinetic equation.
3.2 Macroscopic prediction and microscopic correction
- Macroscopic W updates hold the particle flux fixed, while microscopic P updates hold the wave fixed; accepted particle corrections synchronize conservative moments and refresh the decomposition.The particle solve keeps the collision remainder explicit and bounds relaxation stiffness through its implicit rate.Source: §3 · p. 12, §3.2 · p. 18
3.1 Discrete formulation and residuals
- The wave uses a cell-local characteristic horizon: it adapts from particle-like transport on resolved mean-free-path meshes to Navier–Stokes transport on larger cells, without a switching criterion.The computations use Cℓ = 1 across test cases.Source: §3.1 · pp. 13–14
- Compatible velocity-resolved and macroscopic fluxes, together with conservative collision moments, link the particle and macroscopic residuals so that a vanishing particle residual implies a vanishing macroscopic residual.The two monitored residual norms capture both conservative and nonequilibrium residual components.Source: §3.1 · pp. 15, 17
- The closed discrete steady system RU = 0, RP = 0 is the discrete kinetic equation when moment compatibility is supplied by acceptance.The particle retains the complement omitted by the wave, while the horizon and averaging interval affect transport discretization rather than iteration controls.Source: §3.1 · p. 16
3.2 Macroscopic prediction and microscopic correction
- The endpoint correction evaluates the wave endpoint at the macroscopic prediction and adds it to the first particle update, avoiding relaxation toward the old endpoint; later particle updates need no correction.The correction vanishes at a stationary state.Source: §3.1 · p. 17
3.3.1 Macroscopic point relaxation; 3.3.2 Microscopic LU–SGS
- The implementation uses forward–backward point-relaxation pairs for macroscopic updates and node-independent kinetic LU–SGS sweeps for microscopic updates.The kinetic sweeps operate independently across velocity nodes; the macroscopic sweeps do not depend on the velocity grid.Source: §3.3.1 · p. 19, §3.3.2 · p. 20
4 Consistency and convergence analysis
The analysis establishes exact stationary-state consistency under discrete moment compatibility, while its fluid-limit and prediction-transfer arguments are formal for smooth states. Linear model problems yield convergence factors that are compared with numerical histories.
4 Consistency and convergence analysis
- The stationary-state result is exact under discrete moment compatibility, whereas the fluid-limit and prediction-transfer arguments are formal analyses for smooth states.Convergence factors from linear model problems are compared with numerical histories.Source: §4 · p. 20
- The formal fluid-limit and prediction-transfer analyses are restricted to smooth states.Source: §4 · p. 20
- Linear model problems provide convergence factors that are compared with numerical histories.Source: §4 · p. 20
Discrete collision conditions
The analysis establishes exact stationary-state consistency under discrete moment compatibility, a formal smooth fixed-mesh fluid limit, and model-based convergence mechanisms. Its conclusions are qualified by the discrete collision conditions and by the assumptions of the fluid-limit and convergence analyses.
Discrete collision conditions
- The discrete analysis requires conservative moment matching, a Maxwellian fixed point with zero discrete heat flux, the stated Shakhov heat-flux identity, and horizons bounded below independently of ε.The moment and equilibrium identities are enforced or satisfied to velocity-quadrature accuracy; the exact statements assume those identities.Source: Discrete collision conditions · p. 21
4.1 Stationary-state consistency
- A consistent state is stationary exactly when its particle residual vanishes, making the reconstructed distribution a solution of the discrete kinetic equation independent of iteration controls.A state that passes the stopping test is therefore an approximate solution of that equation, whatever the controls.Source: §4.1 · p. 22
4.2 Hydrodynamic consistency and prediction transfer
- On a fixed mesh, the formal smooth ε → 0 limit recovers finite-volume Euler equations at leading order and the gas-kinetic Navier–Stokes scheme with the specified Prandtl number at first order.The first-order reduction is accurate up to terms of order ε∆i.Source: §4.2 · p. 23
- With smooth near-equilibrium states, locally fixed η, bτ = τ, and the stated solve and acceptance conditions, the accepted distribution and conservative variables match the macroscopic prediction to O(ε).An inexact particle solve requires an O(ε) error in its increment for this transfer result.Source: §4.2 · pp. 23–24
4.3 Model convergence analysis
- For resolved long waves on a fixed mesh, CIS needs O(ε^-2) iterations, while the WPD-based iteration contracts by O(ε²) + O(β(η)) per outer iteration, with β decaying exponentially in η.The model bound requires its stated range; when η < ηc ≈ 1.68, a converged W iteration amplifies long waves.Source: §4.3 · p. 25
- The model analysis also treats finite macroscopic updates and an explicit collision remainder; for Shakhov, the heat-flux mode contracts only at its own rate, 1−Pr = 1/3 per update when Pr = 2/3.This convergence floor belongs to explicit treatment of the remainder, which CIS shares, rather than to the decomposition.Source: §4.3 · pp. 26–27
5 Numerical verification
The numerical assessment tests kinetic fidelity, fluid behavior, convergence mechanisms, and efficiency across rarefied and continuum cases. The setup uses common discretization and computing conditions for the principal IWPD–CIS comparisons and defines separate reference and efficiency measures.
5 Numerical verification
- The assessment spans a normal shock, internal flows, and hypersonic external flows, using profiles and distributions for fidelity and convergence histories, model comparisons, and parameter scans for the implicit coupling.Cavity calculations on grids from a published GSIS study provide an additional solver comparison.Source: §5 · p. 27
5.1 Setup and comparison criteria
- The cases cover ε from 1 to 10^-4 and use the reduced Shakhov model by default, with IWPD and CIS on the same mesh, velocity grid, and CFL numbers.Computations use 64 OpenMP threads on one node; the stated residual tolerance is 10^-7 except for the normal shock and specified cases.Source: §5.1 · pp. 27–28
- Reference agreement is measured by relative difference EL2(q), while efficiency uses iteration and wall-time ratios at the comparison residual threshold.Continuum efficiency is also compared using tWP/tNS; the NS reference shares the mesh, reconstruction, macroscopic solver, and CFL numbers but carries no velocity-space unknowns.Source: §5.1 · p. 28
5.2 Kinetic and hydrodynamic behaviour
Across shock, internal-flow, and hypersonic cases, IWPD retains kinetic profiles while recovering continuum behavior on meshes that do not resolve the mean free path. The tests also show faster convergence than CIS in the continuum regime, with NS faster than IWPD for the reported X38 cylinder comparison.
5.2.1 Normal shock
- For the Mach-10 shock at ε = 1, IWPD and CIS normalized density, velocity, and temperature differ by 3.44 × 10−4 to 4.37×10−4; heat flux and normal stress differ by 1.95×10−3.Reduced distributions differ by 6.39×10−4 to 1.66×10−3, while IWPD captures the upstream hot population and downstream shoulder.Source: §5.2.1 · p. 29, Figure 3 · p. 29
5.2.2 Couette flow
- At ε = 10−4, the Couette velocity profile is the linear continuum solution, and IWPD and CIS differ by 8.73 × 10−4 in velocity and 1.56 × 10−5 in temperature.For ε ≥10−2, velocity differences range from 2.96 × 10−8 to 2.94×10−6 and temperature differences from 7.59×10−11 to 9.27×10−9.Source: §5.2.2 · p. 30, Figure 4 · p. 30
5.2.3 Lid-driven cavity
- At Re = 100 and 1000, IWPD speed fields differ from NS by 8.08×10−3 and 1.54×10−2; its centreline profiles differ from Ghia et al. by 1.87 × 10−3 to 1.41 × 10−2, while CIS departs from the benchmark.The meshes do not resolve the mean free path, and the two methods converge to different discrete solutions.Source: §5.2.3 · p. 30, Figure 5 · p. 31
5.2.4 Hypersonic external flows
- At ε = 1, the cylinder’s rarefied shock layer extends several radii upstream, and IWPD and CIS agree on the stagnation-line and wall temperature and pressure profiles, including the maximum pressure of 16.16.IWPD also reproduces the free-stream beam and reflected-wall populations ahead of the cylinder and the two wake lobes, with relative differences of 3.72 × 10−6.Source: §5.2.4 · pp. 32–33, Figure 6 · p. 32, Figure 7 · p. 33
- At ε = 10−4, IWPD and NS stagnation-line and wall profiles coincide across the bow shock and along the body, although the mesh does not resolve the boundary layer.Wall-profile values are taken from cells adjacent to the wall and therefore differ from wall temperature.Source: §5.2.4 · p. 33, Figure 8 · p. 34
- At ε = 10−4, IWPD reaches the residual tolerance in 20 iterations versus 856 940 for CIS, while NS takes 46.3 s versus 121.8 s for IWPD (tWP/tNS = 2.63).IWPD transports the particle on 642 velocity nodes.Source: §5.2.4 · pp. 33–34, Figure 9 · p. 35
- At ε = 10−2, IWPD–CIS distribution differences at two points near the X38 section are 1.69 × 10−4 and 1.82 × 10−4.The X38 convergence follows the cylinder trend, with the gain over CIS growing as collisions become dominant.Source: §5.2.4 · p. 34, Figure 8 · p. 34
5.3 Convergence and computational cost
CIS loses contraction toward the continuum limit, while IWPD converges in tens to hundreds of outer iterations; the observed rates reflect incomplete macroscopic solves and explicit collision remainders. Near the continuum limit, IWPD’s iteration reductions outweigh its added macroscopic work in wall time.
5.3 Convergence and computational cost
- CIS loses contraction toward the continuum limit, whereas IWPD converges in tens to hundreds of outer iterations.Table 3 compares measured contraction factors fitted over residuals from 10^-6 to 10^-3 with model predictions.Source: §5.3 · pp. 35–36, Table 3 · p. 36
- In the continuum regime, incomplete macroscopic convergence and the explicit collision remainder limit IWPD’s contraction relative to the ideal bound.The ideal bound assumes a converged W iteration; the internal-flow measurements are consistent with the inexact macroscopic solve limiting the outer iteration.Source: §5.3 · p. 36
- For internal flows, IWPD factors are close to the inexact-solve estimate; hypersonic-flow factors at ε = 10−4 lie within 0.002 of the Shakhov value 1−Pr = 1/3.The internal-flow differences are about 0.01 at Re = 100 and 0.02–0.04 for Couette flow and the Re = 1000 cavity.Source: §5.3 · p. 36
- Removing the Shakhov remainder eliminates the homogeneous heat-flux mode, and BGK residuals then fall by more than an order of magnitude per outer iteration below 10−1.Both iteration count and wall time decrease; at ε = 1 and 10−2, the cylinder computations converge similarly for the two collision models.Source: §5.3 · pp. 36–37
5.3.1 Computational cost
- Near the continuum limit, IWPD’s wall-time gains over CIS reach 13.56–50.07 for internal flows and 313.23–2118.23 for external flows.The gains are smaller than the iteration-count ratios because IWPD performs multiple W updates per outer iteration.Source: §5.3.1 · p. 37
- At ε = 1, IWPD and NS require similar iteration counts, and IWPD’s added macroscopic work yields St = 0.40–0.74.For continuum internal flows, the reduction in outer iterations outweighs added macroscopic work; cavity wall-time ratios below one reflect different work allocation between IWPD and NS.Source: §5.3.1 · p. 37
5.4 Parameter and mesh dependence
Parameter scans show that solver cost depends on the allocation of W and P updates, while mesh refinement changes the gain as the cell-to-mean-free-path ratio shifts. On matched grids, IWPD has fewer iterations than GSIS in the reported continuum cases but more at ε = 0.075.
5.4 Parameter and mesh dependence
- At ε = 10−4, the cylinder scan’s lowest wall time is 86.97 s with W58/P2; for large W counts, additional P updates shorten runtime.The scan uses a 100 × 100 mesh with PR2 and covers W1–W256 and P1, P2, P4 and P8.Source: §5.4 · p. 39, Figure 11 · p. 39
- Cost is relatively flat near the optimal W count but more sensitive to point-relaxation work; Couette PR4/W64 is about 18 times slower than its optimum.The common external-flow baseline is W120/PR2/P1; the cylinder scan finds W58/P2 costs about a third as much as W115/P1.Source: §5.4 · p. 39
5.4.1 Mesh refinement and comparison with other methods
- For the refined Re = 100 cavity, CIS iterations rise from 43888 to 50071 as η falls from 5.7 to 2.4, while IWPD’s outer iterations rise from 48 to 187 and St drops from 10.89 to 4.12.The lower η reflects that the finer mesh begins to resolve the mean free path; the gain diminishes with refinement.Source: §5.4.1 · p. 40, Table 6 · p. 40
- At Re = 1000, IWPD wall time remains two thirds of NS across the three meshes, while its advantage over CIS remains St = 56.50 and 56.15 on the two finest meshes.Both IWPD outer iterations and NS cycles grow with cells per direction because the point-relaxation solver slows under refinement.Source: §5.4.1 · p. 40, Table 6 · p. 40
- On the GSIS study’s grids, IWPD requires 150 and 431 outer iterations at Re = 100 and 1000, versus GSIS’s 234 and 1410; at ε = 0.075, the counts are 116 and 52, respectively.The comparisons use each study’s own stopping criterion, so they establish iteration-count comparisons, not a wall-time comparison.Source: §5.4.1 · p. 40
6 Conclusions
The closed WPD system is equivalent to the kinetic equation under moment compatibility and separates Maxwellian wave transport from the particle contribution. The residual-based implicit solver has control-independent stationary states, and the computations support its convergence mechanisms and near-continuum efficiency gains.
6 Conclusions
- Under moment compatibility, the closed WPD system is equivalent to the kinetic equation, with wave and particle contributions to Euler and Navier–Stokes fluxes identified by a hydrodynamic expansion.The wave accumulates only the Maxwellian and depends on conservative variables.Source: §6 · p. 41
- The implicit iteration uses a discrete endpoint correction and moment synchronization, and its stationary states satisfy the target discrete kinetic equation independently of iteration controls.The analysis gives the formal fluid limit on a smooth fixed mesh and distinguishes convergence effects from the WPD scale, incomplete macroscopic solves and explicit collision remainders.Source: §6 · p. 41
- One- and two-dimensional BGK and Shakhov results support the analyzed mechanisms, with substantial iteration-count and wall-time reductions near the continuum limit relative to the conventional iterative scheme.Collision-model comparisons, parameter scans and mesh studies show that gains depend on the collision remainder and allocation of macroscopic work.Source: §6 · p. 41
Supplementary information
Online Resource 1 contains additional profiles, velocity distributions, residual histories, parameter studies and Apollo-section results.
Supplementary information
- Online Resource 1 provides supplementary profiles, velocity distributions, residual histories and parameter studies, including Apollo-section results.Supplementary figures, tables and sections use the prefix S.Source: Supplementary information · p. 42
Declarations
The authors report their contributions, funding, disclosures, and research-material availability. They also disclose the use of AI-assisted tools for manuscript editing.
Author contributions
- Chang Liu led methodology, analysis, software, investigation, visualization and the original draft; Kun Xu contributed to methodology and supervision.Both authors contributed to conceptualization, and both participated in funding acquisition and manuscript review.Source: Declarations · p. 42
Competing interests
- The authors declare no known competing financial interests or personal relationships that could have influenced the reported work.Source: Declarations · p. 42
Funding
- The authors report partial support from Chinese national, Beijing, and Hong Kong funding programs.The listed awards include the National Natural Science Foundation of China, the National Key R&D Program of China, and the Hong Kong Research Grants Council.Source: Declarations · p. 42
Other declarations
- Claude and Codex assisted with manuscript organization and language editing; the authors state they reviewed revisions and remain responsible for the article.Data and code are available from the corresponding author on request; ethics approval and consent are not applicable.Source: Declarations · p. 42
Supplementary material: An implicit multiscale solver for steady kinetic equations based on wave–particle decomposition
The supplement compares IWPD with conventional iteration and, in continuum cases, Navier–Stokes solutions across shock, internal-flow, and hypersonic cases. It reports strong convergence gains for the Apollo case and close kinetic-distribution agreement across sampled X38 and Apollo states.
S1 Normal shock
- The Mach-10 normal-shock figure compares reduced distributions at four positions across the shock, from upstream to downstream, for IWPD and CIS.Source: Figure S1 · p. 45
S3 Lid-driven cavity
- The cavity figures compare IWPD fields and centreline profiles with CIS, and show residual histories and control scans across rarefied and continuum settings.At Re = 100, the figure also compares centreline profiles with Ghia et al. benchmark data [1].Source: Figure S7 · p. 48, Figure S6 · p. 48, Figure S8 · p. 49
S4 Circular cylinder
- Cylinder figures compare IWPD with CIS at ε = 1 and 10−2, and with NS at ε = 10−4, including fields, surface quantities, distributions, and convergence.Source: Figure S10 · p. 51, Figure S11 · p. 51, Figure S12 · p. 52
S5 X38 section
- X38 comparisons show macroscopic fields at Ma = 5 for ε = 1 and 10−2, stagnation-line temperature and pressure at ε = 1, centreline reduced distributions at ε = 1, and NS stagnation-line and wall profiles at ε = 10−4.Source: Figure S17 · p. 55, Figure S18 · p. 55, Figure S19 · p. 56
S6 Additional verification: Apollo section
- Across sampled X38 and Apollo distributions at ε = 1, IWPD–CIS differences ranged from 2.89 × 10−6 to 3.23 × 10−5.Apollo profiles show an asymmetric shock layer at the stated angle of attack, while rarefied distributions include incident and wall-reflected populations and two wake lobes.Source: S6 Additional verification: Apollo… · p. 58
- At ε = 10−4, IWPD reached residual 10−6 in 19 iterations versus 573 864 for CIS; IWPD and NS took 35.47 and 8.89 s, respectively.The corresponding BGK run reached tolerance in 8 iterations.Source: S6 Additional verification: Apollo… · p. 58, Figure S30 · p. 63