1D Euler equations (fluid dynamics) with NDSolve

This is a classical shock-tube problem in which a initially diaphragm separates a hi-pressure, high-density region from one of lower pressure and density. The classical exact solution has multiple discontinuities, a shock wave and a contact-surface (density discontinuity) that propagate to the right, and a continuous rarefaction wave traveling into the high-pressure gas.

The discontinuities cannot be captured by a pseudo spectral scheme,which produces the disastrous oscillations. This also occurs with many other spatial differencing schemes that have too little inherent dissipation. See the book "Difference Methods for Initial Value Problems" by Richtmeyer and Morton (1967), which discusses methods for adding "artificial viscosity" to the Euler equations. This enables the discontinuities to be successfully smeared out over multiple grid points. There are many more modern methods that accomplish the same purpose. Try Googling the terms "Roe Scheme, ENO Scheme, Flux-Difference Splitting."

In short, Mathematica's NDSolve[] command is not yet able to handle directly the Euler equation flows that have discontinuous solutions, and this includes all mixed subsonic/supersonic flows such as your case. Theoretically, it should be able to handle subsonic flows, which are always continuous. Unfortunately, I have found to my dismay that NDSolve[] is not yet sophisticated enough to permit one to use the type of characteristics-based inflow and outflow boundary conditions that are an absolute necessity for generating stable, accurate subsonic flow solutions. I posted a question on this site two days ago, but the question was Closed Out (removed) by other more senior users who deemed it inappropriate. That question was titled "NDSolve rejects a common non-reflective subsonic outflow boundary condition for the 1-D Euler equations".

Update

I've coded a general-purpose function hypersolver to resolve this type of problem.

Here's a solution by adding a very small viscosity μ to Euler equations (notice Euler equations can be seen as particular Navier–Stokes equations with zero viscosity and zero thermal conductivity):

γ = 1.4; p[x_, t_] = (γ - 1) (Ε[x, t] - 1/2 ρ[x, t] v[x, t]^2); eqs = With[{μ = 5 10^-4, ρ = ρ[x, t], v = v[x, t], p = p[x, t], Ε = Ε[x, t]}, {D[ρ, t] + D[ρ v, x] == 0, D[ρ v, t] + D[ρ v^2 + p, x] == μ D[v, x, x](* <- The key point is here *), D[Ε, t] + D[v (Ε + p), x] == 0}]; ρ0[x_] = Piecewise[{{0.125, 0.5 < x <= 1}}, 1]; v0[x_] = 0; p0[x_] = Piecewise[{{0.1, 0.5 < x <= 1}}, 1]; mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, "MinPoints" -> n, "DifferenceOrder" -> o}} {solρ, solv, solΕ} = NDSolveValue[{eqs, {ρ[x, 0] == ρ0[x], ρ[0, t] == ρ0[0], ρ[1, t] == ρ0[1], v[x, 0] == v0[x], v[0, t] == v0[0], v[1, t] == v0[1], p[x, 0] == p0[x], p[0, t] == p0[0], p[1, t] == p0[1]}}, {ρ, v, Ε}, {x, 0, 1}, {t, 0, 0.1}, Method -> mol[400, 4]]; 

μ should be small enough but not too small, or NDSolve will fail again. A dense enough spatial grid is also necessary. Though the warnings eerri and eerr generated, the obtained solution is in good agreement with the one given in OP's link:

Manipulate[GraphicsRow[ Plot[#, {x, 0, 1}, PlotRange -> {0, 11/10}] & /@ {solρ[x, t], solv[x, t], p[x, t] /. {ρ -> solρ, Ε -> solΕ, v -> solv}}, ImageSize -> 1000], {t, 0, 0.1}]