Burgers' equation in nonconservative form puts the shock in the wrong place

Burgers' equation is the simplest model of a nonlinear wave: each value of $u$ travels at speed $u$, so faster parts catch up with slower ones and the profile steepens until it breaks into a jump, a shock. Past that point the derivative $u_x$ does not exist and the differential equation says nothing. What still holds is the integral form: the amount of $u$ between two points changes by what flows in at one end minus what flows out at the other, with flux $f(u)=u^2/2$. Applied to a box around a jump from $u_L$ to $u_R$ moving at speed $s$, it gives the Rankine-Hugoniot condition $s(u_L-u_R)=f(u_L)-f(u_R)$, which for Burgers is $s=(u_L+u_R)/2$.

The runs here start from a step, $u=1$ for $x<0$ and $u=u_R$ for $x>0$, on $-1/2\le x\le 3/2$, with the two states held in a ghost cell at each end, and time step $\Delta t=h/2$. The Python reference, posts/refs/burgers-in-nonconservative-form-puts-the-shock-in-the-wrong-place.py, runs both schemes separately (output beside it).

Two schemes

A conservative scheme updates each cell by the difference of the fluxes through its two faces: $u_i\leftarrow u_i-(\Delta t/h)(F_{i+1/2}-F_{i-1/2})$. Whatever leaves one cell enters its neighbor, so the total $h\sum u_i$ changes only by the fluxes through the two ends of the domain. Godunov's flux is the flux of the exact solution of the jump between the two neighbors; for Burgers it is $\max (f(\max (u_L,0)),f(\min (u_R,0)))$.

The nonconservative scheme discretizes $u_t+u{u}_x=0$ the obvious way. Every $u$ here is positive, so information comes from the left and upwinding takes the left difference: $u_i\leftarrow u_i-(\Delta t/h)u_i(u_i-u_{i-1})$. On smooth data this is a consistent first-order scheme. At a shock the product $u{u}_x$ is a value that jumps times a derivative that is infinite, and the scheme's answer depends on how it approximates that product.

Where the shock goes

Both runs use 800 cells ($h=0.0025$, 800 steps to $t=1$). The shock position is read from the mass: a step from $u_L$ to $u_R$ at $x_s$ holds $u_L(x_s+1/2)+u_R(3/2-x_s)$, so $x_s$ follows from $h\sum u_i$. The speed is the distance moved between $t=1/2$ and $t=1$, divided by $1/2$, which leaves out the first steps, while the profile forms.

Godunov's shock moves at 0.625 to rounding. The mass it holds changes by exactly the inflow $f(u_L)-f(u_R)=0.46875$ per unit time, and that forces the Rankine-Hugoniot speed, whatever the smeared profile looks like. The nonconservative scheme moves its shock at 0.5132. By $t=1$ it has lost 0.085 of mass that the inflow brought in: about $(u_L-u_R)(0.625-0.5132)$, the mass a step would lose by lagging behind at that rate.

Refining the grid does not help

The last cell repeats the nonconservative run on 400 and 1600 cells, and runs $u_R=0$ as well. The speeds on 400, 800 and 1600 cells agree to $3\times 10^{-6}$ (S_SPREAD): the scheme converges, to a solution with the wrong shock speed. A first-order error would halve with each refinement; this one does not change.

With $u_R=0$ the shock should move at $1/2$, and in the nonconservative scheme its speed is 0. The cell just right of the jump has $u_i=0$, so its update $u_i(u_i-u_{i-1})$ is zero whatever its neighbor holds. It stays 0 at every step, and the jump never advances. Godunov's scheme moves this shock at $1/2$ to rounding.