The Spectral Difference method, hands-on

One polynomial per element

Take a conservation law in one dimension,

and cut the domain into elements. Inside each element the Spectral Difference (SD) method keeps the solution at \(p+1\) solution points and represents it by the polynomial of degree \(p\) that passes through them. The scheme is then of order \(p+1\). A second set of \(p+2\) flux points, the two faces of the element plus \(p\) points in between, is where the flux is evaluated.

One update has four moves, and the figure walks through them for advection, \(f = u\). Drag the circles to change the solution and see how each move responds.

In symbols, with \(f_h\) the flux polynomial of degree \(p+1\), the update at solution point \(x_j\) is

There is no integral to compute and no mass matrix to invert, which is what makes the method cheap per degree of freedom and a good fit for GPUs. The interior flux points sit at the Gauss–Legendre points; with that choice the scheme is linearly stable at every order.

The finite volumes hidden inside

The strip at the bottom of the figure shows the intervals between consecutive flux points. Integrate the update over one of them and the derivative of the flux polynomial integrates exactly to a difference of two flux values:

So an SD scheme is also a finite-volume scheme on \(p+1\) unequal control volumes per element, with \(\bar u_i\) the average of the polynomial over each of them, and with an unusually accurate flux at their faces. This second reading is what makes the repair below possible: a flux at one face can be swapped for a more robust one and the scheme remains conservative.

What high order buys

For smooth solutions the error of a scheme of order \(q\) falls like \(h^{q}\). The test below advects a sine wave once around a periodic domain and measures the error for six orders at the same four numbers of degrees of freedom. At 192 degrees of freedom the first-order scheme still has an error of 2 × 10⁻²; the eighth-order one has reached 10⁻¹³, the limit of double precision.

What goes wrong at a discontinuity, and the repair

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 usual cure limits the solution everywhere before it is used. We do the opposite and limit a posteriori:

1. Compute a candidate update with the high-order fluxes, as if nothing could go wrong.
2. Check every control volume of the candidate. Its new value must stay within the range of its neighbours' old values, unless it sits at a smooth extremum, which is detected separately and left alone. For the Euler equations, density and pressure must also remain positive.
3. For the control volumes that fail, replace the fluxes at their faces by those of a robust second-order finite-volume scheme, and redo the update there.

Only the cells that need help are touched, so the solution keeps its high order wherever the flow is smooth. The fallback explainer takes these three steps apart, one control volume at a time.

Things to try

• The overshoot. Press Play with limiting off. The Gaussian is carried almost exactly, and the square wave grows ripples at both edges that persist.
• The repair. Switch limiting on. The ripples disappear. The red marks show which control volumes were recomputed: the ones next to the jumps and no others. Add elements and their number stays about the same while the total grows. The Gaussian keeps its height because its peak is recognised as a smooth extremum.
• Same cost, different order. Choose the sine profile. Order 2 with 32 elements and order 8 with 8 elements both use 64 degrees of freedom. After one period their errors differ by four orders of magnitude.
• Order 1. At order 1 the scheme is the first-order upwind method. Everything is smeared, and nothing ever overshoots.

The demo solves linear advection with a third-order strong-stability-preserving Runge–Kutta integrator and applies the check at every stage. The detection (a relaxed discrete maximum principle with smooth-extremum detection) and the fallback (MUSCL with a minmod slope on the control volumes) follow the reference implementation. The production codes add a blending of the two fluxes around flagged cells and positivity checks on density and pressure.

From this page to the papers

The same construction carries over to systems of equations and to two and three dimensions, where the control volumes become a tensor product of the one-dimensional ones.

• Paper I applies it to the Euler equations in one and two dimensions, with shocks, contact discontinuities and instabilities. Velasco-Romero, Han Veiga & Teyssier (2023)
• Paper II takes it to low Mach numbers, where most shock-capturing schemes dissipate the flow they are meant to follow, with a low-Mach Riemann solver and a well-balanced formulation for stratified atmospheres. Velasco-Romero & Teyssier (2025)
• Paper III adds viscosity and diffusion at arbitrary order. Velasco-Romero & Teyssier (2026)

Where the limiter acts in this run, and in a strong-shock problem, is shown in the fallback explainer: it never touches more than 0.5% of the cells.

The codes are public: spd in Python, which runs on CPUs and GPUs, and spd_K in C++ with Kokkos. There is also a set of notebooks that build the method from scratch.