Centralised Gradient-Based Reconstruction (C-GBR)

• C-GBR is a high-order finite-difference method that uses cell-center values and gradient-based polynomial reconstructions to achieve spectral-like accuracy in smooth regions.
• It centralizes non-acoustic wave reconstructions by averaging left- and right-biased interpolations while applying selective upwinding for acoustic modes to minimize numerical dissipation.
• Coupled with advanced shock sensors and the HLLC Riemann solver, C-GBR maintains robust, grid-independent performance in turbulent and hypersonic flow regimes.

Centralised Gradient-Based Reconstruction (C-GBR) schemes are a class of high-order finite-difference methods for the numerical simulation of compressible flows. These schemes employ compact stencils and gradient-based polynomial reconstruction to achieve spectral-like resolution in smooth regions, while using nonlinear, monotonicity-preserving limiting and shock sensors to maintain stability and sharpness near discontinuities. C-GBR methods are formulated in characteristic space and are tightly integrated with modern approximate Riemann solvers, notably the HLLC flux splitting algorithm. Designed to minimize numerical dissipation away from shocks, C-GBR methods have demonstrated unconditional robustness in challenging regimes, including turbulent boundary layers and hypersonic shock–boundary layer interactions, while remaining competitive in computational performance with artificial-diffusivity-based schemes on modern hardware (Kumar et al., 22 Nov 2025, Chamarthi, 2023, Hoffmann et al., 18 May 2024).

1. Mathematical Formulation and Reconstruction Strategy

C-GBR reconstructs high-order interface states at each grid cell by combining cell-centre point values and high-order finite-difference gradients. The method constructs a two-moment (typically up to second derivative) Legendre polynomial expansion for each variable $\Phi$ in cell $j$:

$$\Phi(x) = \Phi_j + \frac{x-x_j}{\Delta x}\,\Phi'_j + \frac{3\kappa}{2\Delta x^2}\,\Phi''_j \left((x-x_j)^2 - \frac{\Delta x^2}{12}\right),\quad \kappa = \frac{1}{3}$$

Stencils typically use nine points for seventh-order accuracy, or seven/five points for lower-order variants (e.g. MEG6, MIG4). Derivatives are computed via explicit, compact, or implicit centered difference formulas. Crucially, all reconstructions assume that stored cell-center values are pointwise (not cell-averaged), enabling high-order accuracy with compact stencils (Chamarthi, 2023, Kumar et al., 22 Nov 2025).

Centralization is a key innovation: interface states for all but acoustic waves are computed by averaging left- and right-biased interpolations, eliminating the numerical viscosity associated with upwinding for these wave families. For the $b$th characteristic field:

$$C^{C}_{i+1/2,b} = \frac{1}{2} \left( C^{L,\rm Linear}_{i+1/2,b} + C^{R,\rm Linear}_{i+1/2,b} \right)$$

For acoustic waves ($b=1,5$), upwinding and nonlinear limiting remain essential for stability. Entropy/contact and shear/vortical waves ($b=2,3,4$) use centralized interpolation, with monotonicity-preserving limiting activated in regions flagged by shock sensors (Hoffmann et al., 18 May 2024, Kumar et al., 22 Nov 2025).

2. Characteristic-Based Limiting and Discontinuity Detection

Reconstruction is performed in local characteristic space, defined by the right and left eigenvectors of the Jacobian of the normal-to-face convective flux ($\mathbf{A}_n = \mathbf{R}\Lambda\mathbf{L}$). Both primitive and conservative variable transforms have been implemented, with conservative-variable transformations (MEG-C-CONS) systematically providing greater accuracy in nonlinear tests (Hoffmann et al., 18 May 2024).

CRITICAL LIMITERS:

• Monotonicity-Preserving (MP) Limiter: For contact/entropy waves, the limiter of Suresh & Huynh is used, enforcing:

$$\left( C^{L,\rm Linear}_{i+1/2,b} - C_{i,b} \right) \cdot \left( C^{L,\rm Linear}_{i+1/2,b} - C_{i+1,b} \right) \leq \epsilon_{\rm MP}$$

with $\epsilon_{\rm MP}\approx 10^{-20}$ to detect new extrema.

• Ducros Pressure–Vorticity Sensor: For shock detection in vortical (shear/acoustic) waves:

$$\Omega^d = \frac{ |\nabla \cdot u|^2 }{ |\nabla \cdot u|^2 + |\nabla \times u|^2 + \varepsilon }$$

with $\Omega^d > 0.01$ (typical) indicating a strong shock. MP interpolation is then applied locally (Kumar et al., 22 Nov 2025).

3. Integration with HLLC Approximate Riemann Solver

C-GBR methods are coupled to the HLLC approximate Riemann solver for inviscid flux computation. After reconstructing half-cell interface states, the following steps are executed:

1. Transformation back to conservative variables $\mathbf{Q}^{L,R}$.
2. Computation of directional velocities $U^{L,R}$ and sound speeds $c^{L,R}$, followed by Roe-like averages.
3. Determination of eigenvalue bounds $S^L$, $S^R$, and intermediate speed $S_*$ via jump conditions:

$$S_* = \frac{ \rho^R U^R(S^R - U^R) - \rho^L U^L(S^L - U^L) + (p^L-p^R)\|\widetilde{\boldsymbol{\xi}}\|^2 } { \rho^R(S^R - U^R) - \rho^L(S^L - U^L) }$$

1. Computation of star-region conservative states and corresponding fluxes as specified by HLLC schemes (Kumar et al., 22 Nov 2025).
2. Viscous fluxes are computed by standard high-order central schemes, typically reusing gradients already calculated for inviscid reconstruction (Chamarthi, 2023, Hoffmann et al., 18 May 2024).

4. Algorithmic Workflow and Computational Aspects

Each Runge–Kutta stage of the C-GBR solver consists of:

1. High-order gradient extraction for primitive or conservative variables.
2. Transformation to characteristic space and high-order half-cell reconstructions.
3. Centralization and monotonicity-preserving limiting according to wave family and local shock sensors.
4. HLLC Riemann solver flux evaluation, followed by standard central finite-difference viscous fluxes.
5. TVD third-order Runge–Kutta advancement and enforcement of boundary conditions (Kumar et al., 22 Nov 2025).

Numerical efficiency is competitive. On NVIDIA A100, the explicit sixth-order MEG6 variant requires 14–17.4 ms per RK stage for problems of $9$–$17$ million cells, compared to $7.5$–$14.0$ ms for LAD–E4/C6 (artificial-diffusivity schemes). LAD–E4 is $1.9$–$2.3\times$ faster; LAD–C6 is $1.17$–$1.43\times$ faster than MEG6/MIG4, while GPU acceleration achieves $\sim2\times$ speedup over CPU (Kumar et al., 22 Nov 2025).

5. Accuracy, Robustness, and Performance Benchmarks

C-GBR consistently delivers high-resolution solutions in both smooth and discontinuous regimes:

• 1D shock tubes (Sod, Lax, Le Blanc, Shu–Osher): C-GBR resolves shocks in 2–3 cells with negligible overshoot; LAD schemes show spurious oscillations, failing in strong rarefactions (Le Blanc) (Kumar et al., 22 Nov 2025, Chamarthi, 2023).
• 2D and 3D inviscid/viscous tests: C-GBR resolves fine-scale rollups in Richtmyer–Meshkov instability and shock–vortex problems, showing less dissipation than weighted compact nonlinear schemes (WCNS7M) as confirmed by Fourier spectral analysis (Chamarthi, 2023, Hoffmann et al., 18 May 2024).
• 3D turbulent flows (LES/WMLES/TBL, hypersonic ramps, SBLIs): C-GBR (especially MEG-C-CONS) captures separation, reattachment, turbulence onset, and wall heating in agreement with DNS/experiment, even for Mach 7.2 ramp flows. Simulations remain stable for $\mathrm{CFL} \le 0.2$ with no ad hoc filtering. LAD diverges for strong shock–separation (Kumar et al., 22 Nov 2025, Hoffmann et al., 18 May 2024).

In all tests, C-GBR maintains order of accuracy (e.g., seventh order in isentropic vortex, fourth order in pure advection), preserves small-scale enstrophy and retains low dissipation, especially in conservative-variable characteristic space.

6. Hybrid Strategies and Practical Recommendations

For flows dominated by smooth turbulence and moderate shocks (e.g., up to Mach $\sim$3), pure LAD suffices and yields a $\sim1.2\times$ speedup. For strong shock–boundary layer interactions, including compression ramps and triple-points, C-GBR is recommended for its robustness and ability to prevent divergence.

Hybrid solvers activate C-GBR selectively: LAD is applied by default, but when local density or pressure falls below threshold ($\rho<0.5\rho_{\mathrm{ref}}$ or $p<0.5p_\infty$), C-GBR is locally enabled. This approach achieves up to $1.67\times$ speedup versus global MEG6 without loss of stability (Kumar et al., 22 Nov 2025). To maximize computational efficiency, full characteristic limiting is restricted to “troubled” cells as identified by the wave-appropriate sensor.

7. Significance and Impact in High-Speed Compressible Flows

C-GBR methods—by centralizing all non-acoustic wave reconstructions and restricting upwinding (and thus dissipation) to acoustic modes—achieve a combination of shock-capturing capability, spectral accuracy, and low turbulence damping unattainable in traditional upwind or artificial diffusion approaches. The methods are readily adapted to high-order, body-fitted, curvilinear grids and are compatible with wall-modeled LES strategies for hypersonic boundary layer transition (Hoffmann et al., 18 May 2024).

Recent applications demonstrate that the use of conservative-variable characteristic transforms further enhances accuracy in both smooth and shock-dominated flows, outperforming primitive variable transforms in benchmark studies. The approach is robust, free-stream preserving, and achieves grid-independent resolution of turbulence onset, separation bubbles, and heat transfer in regimes where classical artificial-diffusivity or upwind schemes fail or require severe parameter tuning.

A plausible implication is that C-GBR, particularly in hybrid or selectively limited forms, will become a standard tool for high-fidelity, cost-efficient simulation of complex turbulent flows with strong shock–boundary layer interactions in aerodynamic, aerospace, and astrophysical applications (Kumar et al., 22 Nov 2025, Chamarthi, 2023, Hoffmann et al., 18 May 2024).