--- title: Shock capturing summary: Flows with shocks are stabilised by an artificial viscosity that a smoothness sensor switches on only in the elements containing a shock. The term lives inside the differentiated residual, so Newton convergence and adjoint derivatives stay exact. section: theory order: 10 level: 2 keywords: [shock capturing, artificial viscosity, sensor, transonic, Persson-Peraire] --- 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. **[L2 Practice] Using it** ```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. **[L2 Practice] How it works, in brief** 1. **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. 2. **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. 3. **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](../examples/transonic.md)). 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. **[L3 Mechanics] Sensor** 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. **[L3 Mechanics] Viscosity** 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](viscous-br2.md) as the physical viscous terms: lifted gradients, weak divergence, averaged interface flux. No flux crosses domain boundaries. **[L3 Mechanics] Why it is differentiable, and why that matters** 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](implicit-solver.md)), 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. **Table: Adjoint derivatives through the shock** (Example: 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 | Provenance: gyre 0.4.0; revision `39979fcf90` (commit `39979fcf902ad6116eb550b881b4278d4599c8e2`); 2026-11-12T11:37:50-0500; device gpu · NVIDIA GeForce RTX 3060 Ti; quick suite. **Revision check: Same code as this build.** Produced by 39979fcf902ad6116eb550b881b4278d4599c8e2. Every commit between it and the commit this site was built from (ceb1f4b0bed41de53e586df10858079f5a2bcdd8) changes only documentation text, evidence records or bookkeeping: the code is the same. Full record: [example_transonic.json](../data/verification/example_transonic.json). **[L3 Mechanics] Verified properties** 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](viscous-br2.md#penalty-parameter), 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](../examples/transonic.md) checks the flow and its derivatives. > **Warning: Limits** > > The viscosity is constant within an element and discontinuous between elements, which is > adequate for the weak-to-moderate shocks of transonic aerodynamics. It has been exercised on > one transonic aerofoil case at $p = 2$. Strong shocks, supersonic free streams and higher > degrees are untested; positivity is not enforced.