I have implemented a small compressible Euler solver, discretizing in primitive variables (rho, u, v, p) with standard Galerkin FEM P1 triangular elements, and mixed isotropic and anisotropic/streamline diffusion for stabilization.
I have verified the solver on several test cases with shock tube, double shocks, and flow over bumps, so it seems to work fine except for some oscillations (as the artificial diffusion isn't shock capturing or monoticity preserving).
Now I'm trying to validate a simple 2D static oblique shock test case for a unit square which should be simple and has an analytical solution. With top/left boundary conditions Ma_in = 2, rho_in = 1, u_in = cos(10 deg), v_in = sin(10 deg), p_in = 0.1786, bottom slip (v = 0), the analytic solution should be a shock at 29.3 degrees and Ma = 1.64 at the outlet (right boundary). My solver converges just fine and resolves the shock, however the solution is inverted after the shock, Ma 2.4, density 0.65 instead of 1.46 etc.
I cannot seem to find the error, especially as the dual shock problem solves fine. Could this be to using non-conservative/primitive variables? Any ideas at all would be helpful.