Fixed bug in HLLC Riemann solver flux calculation.

This fixes #216 (closed) (again, but now all cases).

I basically overlooked a tiny detail in the real GIZMO and changed the order of some operations that cannot be interchanged. Square test is now verified to work with all 3 possible Riemann solvers.

src/riemann/riemann_hllc.h

@@ -145,11 +145,13 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
      we add the extra velocity flux due to the absolute motion of the fluid
      similarly, we need to add the energy fluxes due to the absolute motion */
   v2 = vij[0] * vij[0] + vij[1] * vij[1] + vij[2] * vij[2];
+  // order is important: we first use the momentum fluxes to update the energy
+  // flux and then de-boost the momentum fluxes!
+  totflux[4] += vij[0] * totflux[1] + vij[1] * totflux[2] +
+                vij[2] * totflux[3] + 0.5 * v2 * totflux[0];
   totflux[1] += vij[0] * totflux[0];
   totflux[2] += vij[1] * totflux[0];
   totflux[3] += vij[2] * totflux[0];
-  totflux[4] += vij[0] * totflux[1] + vij[1] * totflux[2] +
-                vij[2] * totflux[3] + 0.5 * v2 * totflux[0];
 }
 
 #endif /* SWIFT_RIEMANN_HLLC_H */