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}] 

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}] 

Here's a solution by adding a very small viscosity μ to Euler equations:

γ = 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}] 

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}]