Intermittent Sub-grid Wave Correction from Differentiated Riemann Variables

We introduce a low-cost every-$K$-step correction for one-dimensional Euler computations. The correction uses differentiated Riemann variables (DRVs) – characteristic derivatives that isolate the left acoustic wave, the contact, and the right acoustic wave – to locate the current wave packet, sample the surrounding constant states, perform a short Newton update for the intermediate pressure and contact speed, and conservatively remap a sharpened profile back onto the grid. The ingredients are elementary – filtered centered differences, local state sampling, a single Newton step, and conservative cell averaging – yet the effect on accuracy is disproportionate. On a long-time severe-expansion benchmark ($N=900$, $t=0.4$), intermittent correction drives the intermediate-state errors from $O(10^{-2})$ to $O(10^{-13})$, i.e. to machine precision. On a long-time LeBlanc benchmark ($N=800$, $t=1$), the method crosses a qualitative threshold: one-shot final-time reconstruction fails entirely (shock location error $2.7\times 10^{-1}$), whereas correction every three steps recovers an almost exact sharp solution with contact and shock positions accurate to machine precision. The same detector-and-Newton mechanism handles two-shock and two-rarefaction packets without case-specific logic, with plateau improvements of four to sixteen orders of magnitude. In an unoptimized Python prototype the wall-clock overhead is below a factor of two even on the most aggressively corrected benchmark. To our knowledge, no comparably lightweight fixed-grid add-on has been shown to recover this level of coarse-grid accuracy on the long-time LeBlanc and related near-vacuum problems.

💡 Research Summary

The paper introduces a lightweight, intermittent sub‑grid correction technique for one‑dimensional compressible Euler simulations. The method is built on differentiated Riemann variables (DRVs) – the characteristic derivatives ˙w, ˙s, and ˙z – which separate the left acoustic wave, the contact discontinuity, and the right acoustic wave into distinct spikes when evaluated with filtered centered differences. At every K‑th time step the algorithm proceeds as follows:

1. Convert the current cell‑averaged conserved variables (ρ̄, ρū, Ē) into primitive variables (ρ, u, p) and compute algebraic DRV surrogates.
2. Apply a Gaussian filter to the DRV fields and detect the positions of the rarefaction head/tail, the contact, and any shock by examining the sign and magnitude of the spikes. This yields an approximate 1‑2‑3 wave geometry (Xₙʳʰ, Xₙʳᵗ, Xₙᶜ, Xₙˢ).
3. Sample the left, left‑star, right‑star, and right plateau states from the resolved solution.
4. Form a one‑step Newton closure for the star pressure p* and contact speed u* using the classical pressure‑wave function F(p*) = f_L(p*) + f_R(p*) + u_R – u_L = 0, where f_k(p) are the exact acoustic wave functions for each side. The initial guess p⁽⁰⁾* is the average of the sampled left‑star and right‑star pressures; a single Newton update yields p⁽¹⁾* and u⁽¹⁾*.
5. Reconstruct a sharp, self‑similar profile at the current time tₙ based on the recovered p* and u*. The profile is piecewise constant on the plateaus and follows the analytic rarefaction fan formulas inside each detected fan.
6. Replace the existing cell averages with the averages of this reconstructed profile, thereby feeding sub‑cell information back into the subsequent time steps.

Because the correction is applied intermittently (every K steps), the Newton update is effectively distributed over time, acting as a “distributed Newton refinement.” In practice only a few corrections are needed to drive the star‑state pressure from O(10⁻²) down to machine precision (≈10⁻¹³). The method requires no exact Riemann solver in the main loop; the one‑step Newton iteration suffices.

The authors validate the approach on several demanding benchmarks:

• Long‑time severe expansion (γ=1.4, N=900, t=0.4, K=50). The intermediate velocity and pressure errors collapse from 5.7×10⁻² and 2.4×10⁻⁴ to 2.4×10⁻¹³ and 1.8×10⁻¹⁵, essentially machine precision, with negligible wall‑clock penalty.
• Long‑time LeBlanc problem (γ=5/3, N=800, t=1, K=3). A single final‑time post‑processor fails dramatically (shock location error ≈0.27). Intermittent correction every three steps restores both shock and contact positions to machine precision and reduces the L¹ velocity error by four orders of magnitude.
• Two‑shock collision (1‑S/2‑C/3‑S) and two‑rarefaction expansion (1‑R/2‑C/3‑R) (Toro Test 4). The same code path, without any case‑specific logic, yields plateau‑state improvements of five to sixteen orders of magnitude.
• Additional tests include a non‑interacting double‑sod configuration and an extreme near‑vacuum case (Mach 1.5×10⁵) where the contact and shock occupy less than one‑fifth of a cell, confirming the robustness of the detection mechanism.

Performance-wise, an unoptimized Python prototype incurs a wall‑clock overhead ranging from 0.93× (actually faster because the correction reduces the total number of Euler steps) to 1.84× on the most aggressively corrected benchmark. Thus the method is computationally affordable.

Conceptually, the authors draw an analogy with large‑eddy simulation (LES): LES filters the Navier–Stokes equations and models the effect of unresolved eddies on the resolved scales, whereas the DRV correction filters the resolved Euler solution, infers the hidden wave geometry, and injects a deterministic correction back into the evolution. Unlike LES, the sub‑grid structure here is fully deterministic because the 1‑D Riemann problem has a unique analytic solution once the outer states and wave pattern are known.

In summary, the paper demonstrates that a modest, deterministic sub‑grid correction—based on DRV detection, simple plateau sampling, a single Newton step, and conservative remapping—can dramatically improve the accuracy of standard fixed‑grid shock‑capturing schemes. It achieves near‑exact reconstruction of shock, contact, and rarefaction structures on coarse grids, turning otherwise failed long‑time simulations into almost exact solutions, all with minimal computational cost. This makes the technique a powerful, low‑overhead add‑on for existing high‑order Euler solvers.