Shock capturing

A polynomial cannot represent a discontinuity. Entropy-conservative volume terms and a dissipative interface flux remove the aliasing instability, but at a shock the solution still oscillates, and nothing in the scheme keeps density and pressure positive. Gyre adds a small amount of physical-looking diffusion exactly where the solution stops being resolved, and nowhere else.

bash

.venv/bin/gyre run naca0012_transonic            # preset: M = 0.8, α = 1.25°
.venv/bin/gyre run naca0012_euler --set flow.mach=0.8 --set discretization.shock_capturing=true

discretization.shock_capturing enables the term; discretization.av_coef scales it (default 1). Smooth flows are unaffected: the sensor is exactly zero for a uniform state and far below threshold for resolved solutions.

Detect. In each element the density is expanded in Legendre modes. If a noticeable part of its energy sits in the highest modes, the element is not resolving the solution — the signature of a shock. This gives one smooth sensor value per element.
Diffuse. Where the sensor exceeds a threshold, a Laplacian artificial viscosity is switched on, scaled with the element size divided by the polynomial degree and with the local wave speed, and with av_coef. It spreads the shock over a few nodes inside the element instead of letting the polynomial oscillate.
Differentiate. Sensor and viscosity are smooth functions of the state and part of the residual, so the Newton linearisation and the adjoint include them exactly. This is why a transonic case converges to machine precision and why its gradients agree with finite differences (transonic example).

What it is not: a limiter, a positivity guarantee or an entropy-stable shock treatment. It is verified on a transonic aerofoil with weak shocks and untested for strong ones; the limits are listed at the end of this page.

In every element the density is transformed to an orthonormal Legendre basis. The sensor is the fraction of its energy carried by the highest polynomial shell (all modes with at least one index equal to $p$ ):

S_e = \frac{\sum_{\text{modes with } \max(i,j,k) = p} \hat\rho^{\,2}}{\sum_{\text{all modes}} \hat\rho^{\,2}}, \qquad s_e = \log_{10} S_e .

For a smooth, resolved solution the spectrum decays rapidly and $s_e$ is very negative; across a shock the spectrum is flat and $s_e$ approaches zero. The numerator is evaluated as a direct sum of squares — computing it as “total minus low modes” lets round-off make it negative in near-uniform far-field elements, and the logarithm of a negative number poisons the whole linearisation.

The element-wise coefficient ramps smoothly from zero to its maximum,

\varepsilon_e = \varepsilon_{\max}\cdot \begin{cases} 0 & s_e \le s_0 - \kappa \\ \tfrac12\left(1 + \sin\dfrac{\pi (s_e - s_0)}{2\kappa}\right) & |s_e - s_0| < \kappa \\ 1 & s_e \ge s_0 + \kappa \end{cases} \qquad \varepsilon_{\max} = c_\text{av}\,\frac{h_e}{p}\,\overline{(|\mathbf u| + c)}_e ,

with threshold $s_0 = -(1.5 + 4\log_{10} p)$ , half-width $\kappa = 1$ , $h_e$ the element size and the bar an element mean. The element size is the square root of the element area in 2-D and the cube root of the volume in 3-D, both taken from the nodal mean of the Jacobian. On a spanwise extrusion of a 2-D section it is the square root of the in-plane area, so that an extruded case carries, for the same state, exactly the viscosity of its section, independent of span and nz. The scaling makes the shock about one element wide. The diffusive flux $\varepsilon_e \nabla\mathbf U$ acts on all conserved variables and is discretized with the same BR2 machinery as the physical viscous terms: lifted gradients, weak divergence, averaged interface flux. No flux crosses domain boundaries.

The sensor, the ramp (which is $C^1$ ) and the viscosity are ordinary JAX code inside the residual. Therefore:

the Newton linearisation includes the dependence of $\varepsilon_e$ on the solution — the solver reduces the residual of a transonic case by ten orders, at the linear rate of its inexact Newton iteration (implicit solver), rather than stalling in the limit cycle that frozen or non-smooth limiters are known for;
the element blocks of the preconditioner remain exact (the viscosity of an element depends only on its own state; the neighbour’s enters through the neighbour-role argument);
the discrete adjoint contains the shock-capturing term, so total derivatives are those of the function that is actually computed.

Adjoint derivatives through the shockExample: transonic aerofoil with shock capturing

derivative	adjoint	finite difference	relative difference
dcl/dalpha	17.6219	17.6219	1.4558e-06
dcd/dalpha	1.31459	1.31459	9.9160e-07
dcl/dmach	4.46992	4.46988	9.9365e-06
dcd/dmach	0.688727	0.688725	3.7501e-06

gyre 0.4.0 · 39979fcf90 · 2026-11-12 11:37 · gpu · NVIDIA GeForce RTX 3060 Ti · quick suite · record (JSON)Same code as this build

Unit tests assert that the viscosity is exactly zero for uniform and smooth resolved data and non-zero only in elements at a steep profile; that free-stream preservation and conservation hold with the term active (2-D and 3-D); that it increases entropy dissipation; that a spanwise extrusion reproduces the viscosity of its 2-D section for different spans, and its residual for a state that is continuous across elements (where the state jumps between elements the two residuals differ through the default penalty, 6 for hexahedra against 4 for quadrilaterals); and that block inverses, JVPs and VJPs remain exact with the sensor partially ramped (tests/test_operator.py). The transonic example checks the flow and its derivatives.