The fallback scheme, hands-on

Why a repair is needed

A polynomial cannot follow a jump. It overshoots on both sides, and the overshoots do not shrink when the mesh is refined. For the Euler equations an overshoot can mean a negative density or pressure, and the simulation stops. The advection demo of the Spectral Difference page shows it at any order above one.

The classical cure acts a priori: before each update, slopes or high modes are reduced wherever a sensor judges the solution too rough. The sensor has to guess where trouble will appear, and what it removes from smooth flow by mistake is lost. The alternative, introduced by Clain, Diot & Loubère (2011) under the name MOOD, acts a posteriori: take the unlimited step, look at what came out, and redo only what is not acceptable.

That needs two ingredients: a robust scheme to fall back on that can be mixed with the high-order one cell by cell, and a test that tells an acceptable value from a bad one.

The finite volumes inside an element

The first ingredient comes for free. A Spectral Difference (SD) element of degree \(p\) evaluates its flux at \(p+2\) flux points, and the \(p+1\) intervals between them are control volumes. Integrating the SD update over one of them gives, exactly,

with \(\bar u_i\) the average of the solution over the control volume and \(f_h\) the flux polynomial (how it is built). So the SD scheme is already a finite-volume scheme on these control volumes, with an unusually accurate flux at their faces. The fallback is a second finite-volume scheme on the same control volumes. Both are differences of fluxes at the same faces, so the flux at any one face can be taken from either scheme and the result is still conservative.

One stage, move by move

Every stage of the time integrator goes through the same five moves. The figure shows them for a square pulse advected to the right, \(f = u\), on eight fourth-order elements.

Things to try

• Follow the pulse. Go to move 3 and press Play. The rejected volumes sit where the smeared edges of the pulse meet the flat parts, and travel with them. A few also appear on the flat parts themselves: there the band has no width, and a ripple of one part in 10⁵ is enough to leave it.
• A smooth profile. Choose the sine. The crest and the trough leave their bands at every step, as they must in order to move. They are marked blue and kept. Nothing is rejected, and the scheme is the plain high-order one.
• A smooth maximum. Choose the Gaussian. Its top is marked blue and kept in the same way. The few red marks are far out in the tails, where the profile is flat to five digits and the band has no width.
• The tolerance. Raise it to 10⁻³: the rejections in the tails of the Gaussian and on the flat parts of the square are gone, and those at the edges of the square become rarer. At 10⁻² and above the edges are let through with a visible overshoot.

The checks

The candidate value of a control volume is rejected if it fails the first test without being excused by the second, or if it fails the third.

• Numerical admissibility. A scheme that only transports the solution cannot create new extrema. With \(\bar u^{\,n}\) the values before the stage and \(\bar u^{\,*}\) the candidate, the test is a discrete maximum principle over the volume and its face neighbours,
  \[ m_i - \varepsilon\,|m_i| \;\le\; \bar u_i^{\,*} \;\le\; M_i + \varepsilon\,|M_i| , \qquad m_i = \min\left(\bar u^{\,n}_{i-1}, \bar u^{\,n}_{i}, \bar u^{\,n}_{i+1}\right), \quad M_i = \max\left(\bar u^{\,n}_{i-1}, \bar u^{\,n}_{i}, \bar u^{\,n}_{i+1}\right) . \]
  The tolerance \(\varepsilon\) keeps round-off and harmless ripples from triggering it; typical values are \(10^{-5}\) to \(10^{-3}\). For the Euler equations the test is applied to density and pressure.
• Smooth extrema. A genuine maximum or minimum of a smooth solution violates that principle whenever it moves, and clipping it would cost the scheme its accuracy there. So a volume that fails the first test is examined again: the change of the candidate's slope towards each neighbour is compared with the change predicted by the local curvature. If the two agree in the volume and in both of its neighbours, the profile is locally a smooth parabola and the value is accepted.
• Physical admissibility. Density and pressure must be positive. A candidate that is not is rejected whatever the other tests say.

The repair

The fallback is a second-order Godunov scheme, MUSCL–Hancock, applied to the control volumes: a slope is built in each one from its neighbours and limited (minmod in the demos), and a Riemann problem is solved at every face. It is diffusive and it does not overshoot, which is what is wanted in the few places where it is used.

At a face of a rejected volume the flux is taken from the fallback; elsewhere it is the SD flux. In general the two are mixed with a weight \(\theta\) defined at the faces,

with \(\theta = 1\) at the faces of a rejected volume. The figure above uses \(\theta = 0\) everywhere else; the production codes let \(\theta\) fall to ¾ and ¼ over the next two faces, so that the scheme does not change abruptly from one volume to the next. The update is then repeated with the mixed fluxes. Because each face still carries a single flux, what leaves one volume enters its neighbour, at every order and with any \(\theta\).

All of this happens at every stage of the time integrator, on the control volumes of every element at once. In two and three dimensions the control volumes are tensor products of the one-dimensional ones and the same moves are applied face by face, direction by direction.

A shock that forms on its own

In a gas, shocks are not put in by hand: they form. The simplest model is Burgers' equation, \(f = u^2/2\), in which each value travels at a speed equal to itself, so the crest of a wave catches up with the trough ahead of it. A sine wave stays smooth until \(t = 1/\pi \approx 0.32\) and carries a shock from then on. A limiter that acts after the fact has nothing to do before that time and a great deal to do after it.

Things to try

• Silence, then the shock. Press Play. The red curve in the strip stays at zero while the wave steepens, and leaves it at \(t = 0.28\), shortly before the shock forms, when the front has become too steep for sixteen elements. From then on only the few control volumes around the shock are recomputed.
• Without the repair. Untick the limiting and play again. Up to the first fallback the two runs are identical; after it the unlimited one rings on both sides of the jump.
• More elements. With 32 or 64 elements the wave stays resolved for longer and the first fallback moves towards \(t = 0.32\). About five volumes are repaired per step whatever the number of elements, so their share halves each time the mesh is doubled.
• Order 2. The smooth-extremum test needs a parabola, so it is not used at order 2. The crest and the trough are then repaired from the first step on, as in any limited second-order scheme.

In real simulations

The same construction runs in the production codes for the Euler and Navier–Stokes equations, and the share of the flow that is repaired stays small. In the two runs below, the control volumes that fell back are marked in red.

The demos solve linear advection and Burgers' equation with a third-order strong-stability-preserving Runge–Kutta integrator and apply the checks at every stage. The detection (the relaxed maximum principle with the smooth-extremum test) and the fallback (MUSCL with a minmod slope on the unequal control volumes) follow the reference implementation. The codes of the papers integrate in time with ADER, correct each of its time slices in the same way, and add the flux blending and the positivity checks described above.

From this page to the papers

• Paper I introduces the scheme for the Euler equations in one and two dimensions. Velasco-Romero, Han Veiga & Teyssier (2023)
• Paper II takes it to low Mach numbers and stratified atmospheres, where the checks are applied to the departure from a hydrostatic equilibrium. Velasco-Romero & Teyssier (2025)
• Paper III adds viscosity and diffusion at arbitrary order. Velasco-Romero & Teyssier (2026)

The codes are public: spd in Python, which runs on CPUs and GPUs, and spd_K in C++ with Kokkos.