diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_viscous_shock.jl
new file mode 100644
index 00000000000..feec3d534d0
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_viscous_shock.jl
@@ -0,0 +1,183 @@
+using OrdinaryDiffEq
+using Trixi
+
+# This is the classic 1D viscous shock wave problem with analytical solution 
+# for a special value of the Prandtl number.
+# The original references are:
+#
+# - R. Becker (1922)
+#   Stoßwelle und Detonation.
+#   [DOI: 10.1007/BF01329605](https://doi.org/10.1007/BF01329605)
+#
+#   English translations:
+#   Impact waves and detonation. Part I.
+#   https://ntrs.nasa.gov/api/citations/19930090862/downloads/19930090862.pdf
+#   Impact waves and detonation. Part II.
+#   https://ntrs.nasa.gov/api/citations/19930090863/downloads/19930090863.pdf
+#
+# - M. Morduchow, P. A. Libby (1949)
+#   On a Complete Solution of the One-Dimensional Flow Equations 
+#   of a Viscous, Head-Conducting, Compressible Gas
+#   [DOI: 10.2514/8.11882](https://doi.org/10.2514/8.11882)
+#
+#
+# The particular problem considered here is described in
+# - L. G. Margolin, J. M. Reisner, P. M. Jordan (2017) 
+#   Entropy in self-similar shock profiles
+#   [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+
+### Fixed parameters ###
+
+# Special value for which nonlinear solver can be omitted
+# Corresponds essentially to fixing the Mach number
+alpha = 0.5
+# We want kappa = cp * mu = mu_bar to ensure constant enthalpy
+prandtl_number() = 3 / 4
+
+### Free choices: ###
+gamma() = 5 / 3
+
+# In Margolin et al., the Navier-Stokes equations are given for an 
+# isotropic stress tensor τ, i.e., ∇ ⋅ τ = μ Δu 
+mu_isotropic() = 0.15
+mu_bar() = mu_isotropic() / (gamma() - 1) # Re-scaled viscosity
+
+rho_0() = 1
+v() = 1 # Shock speed
+
+domain_length = 4.0
+
+### Derived quantities ###
+
+Ma() = 2 / sqrt(3 - gamma()) # Mach number for alpha = 0.5
+c_0() = v() / Ma() # Speed of sound ahead of the shock
+
+# From constant enthalpy condition
+p_0() = c_0()^2 * rho_0() / gamma()
+
+l() = mu_bar() / (rho_0() * v()) * 2 * gamma() / (gamma() + 1) # Appropriate length scale
+
+"""
+    initial_condition_viscous_shock(x, t, equations)
+
+Classic 1D viscous shock wave problem with analytical solution 
+for a special value of the Prandtl number.
+The version implemented here is described in
+- L. G. Margolin, J. M. Reisner, P. M. Jordan (2017)
+  Entropy in self-similar shock profiles
+  [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+"""
+function initial_condition_viscous_shock(x, t, equations)
+    y = x[1] - v() * t # Translated coordinate
+
+    # Coordinate transformation. See eq. (33) in Margolin et al. (2017)
+    chi = 2 * exp(y / (2 * l()))
+
+    w = 1 + 1 / (2 * chi^2) * (1 - sqrt(1 + 2 * chi^2))
+
+    rho = rho_0() / w
+    u = v() * (1 - w)
+    p = p_0() * 1 / w * (1 + (gamma() - 1) / 2 * Ma()^2 * (1 - w^2))
+
+    return prim2cons(SVector(rho, u, 0, p), equations)
+end
+initial_condition = initial_condition_viscous_shock
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+equations = CompressibleEulerEquations2D(gamma())
+
+# Trixi implements the stress tensor in deviatoric form, thus we need to 
+# convert the "isotropic viscosity" to the "deviatoric viscosity"
+mu_deviatoric() = mu_bar() * 3 / 4
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu_deviatoric(),
+                                                          Prandtl = prandtl_number(),
+                                                          gradient_variables = GradientVariablesPrimitive())
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_hlle)
+
+coordinates_min = (-domain_length / 2, -domain_length / 2)
+coordinates_max = (domain_length / 2, domain_length / 2)
+
+trees_per_dimension = (8, 2)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3, initial_refinement_level = 0,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = (false, true))
+
+### Inviscid boundary conditions ###
+
+# Prescribe pure influx based on initial conditions
+function boundary_condition_inflow(u_inner, normal_direction::AbstractVector, x, t,
+                                   surface_flux_function,
+                                   equations::CompressibleEulerEquations2D)
+    u_cons = initial_condition_viscous_shock(x, t, equations)
+    flux = Trixi.flux(u_cons, normal_direction, equations)
+
+    return flux
+end
+
+# Completely free outflow
+function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t,
+                                    surface_flux_function,
+                                    equations::CompressibleEulerEquations2D)
+    # Calculate the boundary flux entirely from the internal solution state
+    flux = Trixi.flux(u_inner, normal_direction, equations)
+
+    return flux
+end
+
+boundary_conditions = Dict(:x_neg => boundary_condition_inflow,
+                           :x_pos => boundary_condition_outflow)
+
+### Viscous boundary conditions ###
+# For the viscous BCs, we use the known analytical solution
+velocity_bc = NoSlip() do x, t, equations_parabolic
+    Trixi.velocity(initial_condition_viscous_shock(x,
+                                                   t,
+                                                   equations_parabolic),
+                   equations_parabolic)
+end
+
+heat_bc = Isothermal() do x, t, equations_parabolic
+    Trixi.temperature(initial_condition_viscous_shock(x,
+                                                      t,
+                                                      equations_parabolic),
+                      equations_parabolic)
+end
+
+boundary_condition_parabolic = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc)
+
+boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_parabolic,
+                                     :x_pos => boundary_condition_parabolic)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions = (boundary_conditions,
+                                                                    boundary_conditions_parabolic))
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+alive_callback = AliveCallback(alive_interval = 10)
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-8
+sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
+            dt = 1e-3, ode_default_options()..., callback = callbacks)
+
+summary_callback() # print the timer summary
diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_viscous_shock.jl
new file mode 100644
index 00000000000..43cdafa3753
--- /dev/null
+++ b/examples/p4est_3d_dgsem/elixir_navierstokes_viscous_shock.jl
@@ -0,0 +1,183 @@
+using OrdinaryDiffEq
+using Trixi
+
+# This is the classic 1D viscous shock wave problem with analytical solution 
+# for a special value of the Prandtl number.
+# The original references are:
+#
+# - R. Becker (1922)
+#   Stoßwelle und Detonation.
+#   [DOI: 10.1007/BF01329605](https://doi.org/10.1007/BF01329605)
+#
+#   English translations:
+#   Impact waves and detonation. Part I.
+#   https://ntrs.nasa.gov/api/citations/19930090862/downloads/19930090862.pdf
+#   Impact waves and detonation. Part II.
+#   https://ntrs.nasa.gov/api/citations/19930090863/downloads/19930090863.pdf
+#
+# - M. Morduchow, P. A. Libby (1949)
+#   On a Complete Solution of the One-Dimensional Flow Equations 
+#   of a Viscous, Head-Conducting, Compressible Gas
+#   [DOI: 10.2514/8.11882](https://doi.org/10.2514/8.11882)
+#
+#
+# The particular problem considered here is described in
+# - L. G. Margolin, J. M. Reisner, P. M. Jordan (2017) 
+#   Entropy in self-similar shock profiles
+#   [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+
+### Fixed parameters ###
+
+# Special value for which nonlinear solver can be omitted
+# Corresponds essentially to fixing the Mach number
+alpha = 0.5
+# We want kappa = cp * mu = mu_bar to ensure constant enthalpy
+prandtl_number() = 3 / 4
+
+### Free choices: ###
+gamma() = 5 / 3
+
+# In Margolin et al., the Navier-Stokes equations are given for an 
+# isotropic stress tensor τ, i.e., ∇ ⋅ τ = μ Δu 
+mu_isotropic() = 0.15
+mu_bar() = mu_isotropic() / (gamma() - 1) # Re-scaled viscosity
+
+rho_0() = 1
+v() = 1 # Shock speed
+
+domain_length = 4.0
+
+### Derived quantities ###
+
+Ma() = 2 / sqrt(3 - gamma()) # Mach number for alpha = 0.5
+c_0() = v() / Ma() # Speed of sound ahead of the shock
+
+# From constant enthalpy condition
+p_0() = c_0()^2 * rho_0() / gamma()
+
+l() = mu_bar() / (rho_0() * v()) * 2 * gamma() / (gamma() + 1) # Appropriate length scale
+
+"""
+    initial_condition_viscous_shock(x, t, equations)
+
+Classic 1D viscous shock wave problem with analytical solution 
+for a special value of the Prandtl number.
+The version implemented here is described in
+- L. G. Margolin, J. M. Reisner, P. M. Jordan (2017)
+  Entropy in self-similar shock profiles
+  [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+"""
+function initial_condition_viscous_shock(x, t, equations)
+    y = x[1] - v() * t # Translated coordinate
+
+    # Coordinate transformation. See eq. (33) in Margolin et al. (2017)
+    chi = 2 * exp(y / (2 * l()))
+
+    w = 1 + 1 / (2 * chi^2) * (1 - sqrt(1 + 2 * chi^2))
+
+    rho = rho_0() / w
+    u = v() * (1 - w)
+    p = p_0() * 1 / w * (1 + (gamma() - 1) / 2 * Ma()^2 * (1 - w^2))
+
+    return prim2cons(SVector(rho, u, 0, 0, p), equations)
+end
+initial_condition = initial_condition_viscous_shock
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+equations = CompressibleEulerEquations3D(gamma())
+
+# Trixi implements the stress tensor in deviatoric form, thus we need to 
+# convert the "isotropic viscosity" to the "deviatoric viscosity"
+mu_deviatoric() = mu_bar() * 3 / 4
+equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu_deviatoric(),
+                                                          Prandtl = prandtl_number(),
+                                                          gradient_variables = GradientVariablesPrimitive())
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_hlle)
+
+coordinates_min = (-domain_length / 2, -domain_length / 2, -domain_length / 2)
+coordinates_max = (domain_length / 2, domain_length / 2, domain_length / 2)
+
+trees_per_dimension = (8, 2, 2)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3, initial_refinement_level = 0,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = (false, true, true))
+
+### Inviscid boundary conditions ###
+
+# Prescribe pure influx based on initial conditions
+function boundary_condition_inflow(u_inner, normal_direction::AbstractVector, x, t,
+                                   surface_flux_function,
+                                   equations::CompressibleEulerEquations3D)
+    u_cons = initial_condition_viscous_shock(x, t, equations)
+    flux = Trixi.flux(u_cons, normal_direction, equations)
+
+    return flux
+end
+
+# Completely free outflow
+function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t,
+                                    surface_flux_function,
+                                    equations::CompressibleEulerEquations3D)
+    # Calculate the boundary flux entirely from the internal solution state
+    flux = Trixi.flux(u_inner, normal_direction, equations)
+
+    return flux
+end
+
+boundary_conditions = Dict(:x_neg => boundary_condition_inflow,
+                           :x_pos => boundary_condition_outflow)
+
+### Viscous boundary conditions ###
+# For the viscous BCs, we use the known analytical solution
+velocity_bc = NoSlip() do x, t, equations_parabolic
+    Trixi.velocity(initial_condition_viscous_shock(x,
+                                                   t,
+                                                   equations_parabolic),
+                   equations_parabolic)
+end
+
+heat_bc = Isothermal() do x, t, equations_parabolic
+    Trixi.temperature(initial_condition_viscous_shock(x,
+                                                      t,
+                                                      equations_parabolic),
+                      equations_parabolic)
+end
+
+boundary_condition_parabolic = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc)
+
+boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_parabolic,
+                                     :x_pos => boundary_condition_parabolic)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions = (boundary_conditions,
+                                                                    boundary_conditions_parabolic))
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+alive_callback = AliveCallback(alive_interval = 10)
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-8
+sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
+            dt = 1e-3, ode_default_options()..., callback = callbacks)
+
+summary_callback() # print the timer summary
diff --git a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
index 3cb4735e132..6b7dd007128 100644
--- a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
+++ b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
@@ -19,7 +19,8 @@ mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
 # Green light that at x = 0 which switches at t = 0 from red to green.
 # To the left there are cars bumper to bumper, to the right there are no cars.
 function initial_condition_greenlight(x, t, equation::TrafficFlowLWREquations1D)
-    scalar = x[1] < 0.0 ? 1.0 : 0.0
+    RealT = eltype(x)
+    scalar = x[1] < 0 ? one(RealT) : zero(RealT)
 
     return SVector(scalar)
 end
@@ -29,7 +30,8 @@ end
 
 # Assume that there are always cars waiting at the left 
 function inflow(x, t, equations::TrafficFlowLWREquations1D)
-    return initial_condition_greenlight(coordinates_min, t, equations)
+    # -1.0 = coordinates_min
+    return initial_condition_greenlight(-1.0, t, equations)
 end
 boundary_condition_inflow = BoundaryConditionDirichlet(inflow)
 
diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl
new file mode 100644
index 00000000000..df6a453ac83
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl
@@ -0,0 +1,176 @@
+using OrdinaryDiffEq
+using Trixi
+
+# This is the classic 1D viscous shock wave problem with analytical solution 
+# for a special value of the Prandtl number.
+# The original references are:
+#
+# - R. Becker (1922)
+#   Stoßwelle und Detonation.
+#   [DOI: 10.1007/BF01329605](https://doi.org/10.1007/BF01329605)
+#
+#   English translations:
+#   Impact waves and detonation. Part I.
+#   https://ntrs.nasa.gov/api/citations/19930090862/downloads/19930090862.pdf
+#   Impact waves and detonation. Part II.
+#   https://ntrs.nasa.gov/api/citations/19930090863/downloads/19930090863.pdf
+#
+# - M. Morduchow, P. A. Libby (1949)
+#   On a Complete Solution of the One-Dimensional Flow Equations 
+#   of a Viscous, Head-Conducting, Compressible Gas
+#   [DOI: 10.2514/8.11882](https://doi.org/10.2514/8.11882)
+#
+#
+# The particular problem considered here is described in
+# - L. G. Margolin, J. M. Reisner, P. M. Jordan (2017) 
+#   Entropy in self-similar shock profiles
+#   [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+
+### Fixed parameters ###
+
+# Special value for which nonlinear solver can be omitted
+# Corresponds essentially to fixing the Mach number
+alpha = 0.5
+# We want kappa = cp * mu = mu_bar to ensure constant enthalpy
+prandtl_number() = 1
+
+### Free choices: ###
+gamma() = 5 / 3
+
+mu() = 0.15
+mu_bar() = mu() / (gamma() - 1) # Re-scaled viscosity
+
+rho_0() = 1
+v() = 1 # Shock speed
+
+domain_length = 4.0
+
+### Derived quantities ###
+
+Ma() = 2 / sqrt(3 - gamma()) # Mach number for alpha = 0.5
+c_0() = v() / Ma() # Speed of sound ahead of the shock
+
+# From constant enthalpy condition
+p_0() = c_0()^2 * rho_0() / gamma()
+
+l() = mu_bar() / (rho_0() * v()) * 2 * gamma() / (gamma() + 1) # Appropriate length scale
+
+"""
+    initial_condition_viscous_shock(x, t, equations)
+
+Classic 1D viscous shock wave problem with analytical solution 
+for a special value of the Prandtl number.
+The version implemented here is described in
+- L. G. Margolin, J. M. Reisner, P. M. Jordan (2017)
+  Entropy in self-similar shock profiles
+  [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+"""
+function initial_condition_viscous_shock(x, t, equations)
+    y = x[1] - v() * t # Translated coordinate
+
+    # Coordinate transformation. See eq. (33) in Margolin et al. (2017)
+    chi = 2 * exp(y / (2 * l()))
+
+    w = 1 + 1 / (2 * chi^2) * (1 - sqrt(1 + 2 * chi^2))
+
+    rho = rho_0() / w
+    u = v() * (1 - w)
+    p = p_0() * 1 / w * (1 + (gamma() - 1) / 2 * Ma()^2 * (1 - w^2))
+
+    return prim2cons(SVector(rho, u, p), equations)
+end
+initial_condition = initial_condition_viscous_shock
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+equations = CompressibleEulerEquations1D(gamma())
+equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu = mu_bar(),
+                                                          Prandtl = prandtl_number(),
+                                                          gradient_variables = GradientVariablesPrimitive())
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_hlle)
+
+coordinates_min = -domain_length / 2
+coordinates_max = domain_length / 2
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                periodicity = false,
+                n_cells_max = 30_000)
+
+### Inviscid boundary conditions ###
+
+# Prescribe pure influx based on initial conditions
+function boundary_condition_inflow(u_inner, orientation::Integer, normal_direction, x, t,
+                                   surface_flux_function,
+                                   equations::CompressibleEulerEquations1D)
+    u_cons = initial_condition_viscous_shock(x, t, equations)
+    flux = Trixi.flux(u_cons, orientation, equations)
+
+    return flux
+end
+
+# Completely free outflow
+function boundary_condition_outflow(u_inner, orientation::Integer, normal_direction, x, t,
+                                    surface_flux_function,
+                                    equations::CompressibleEulerEquations1D)
+    # Calculate the boundary flux entirely from the internal solution state
+    flux = Trixi.flux(u_inner, orientation, equations)
+
+    return flux
+end
+
+boundary_conditions = (; x_neg = boundary_condition_inflow,
+                       x_pos = boundary_condition_outflow)
+
+### Viscous boundary conditions ###
+# For the viscous BCs, we use the known analytical solution
+velocity_bc = NoSlip() do x, t, equations_parabolic
+    Trixi.velocity(initial_condition_viscous_shock(x,
+                                                   t,
+                                                   equations_parabolic),
+                   equations_parabolic)
+end
+
+heat_bc = Isothermal() do x, t, equations_parabolic
+    Trixi.temperature(initial_condition_viscous_shock(x,
+                                                      t,
+                                                      equations_parabolic),
+                      equations_parabolic)
+end
+
+boundary_condition_parabolic = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc)
+
+boundary_conditions_parabolic = (; x_neg = boundary_condition_parabolic,
+                                 x_pos = boundary_condition_parabolic)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions = (boundary_conditions,
+                                                                    boundary_conditions_parabolic))
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+alive_callback = AliveCallback(alive_interval = 100)
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-8
+sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
+            dt = 1e-3, ode_default_options()..., callback = callbacks)
+
+summary_callback() # print the timer summary
diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl
index 024ce2d7e87..fe79a7e0590 100644
--- a/src/equations/compressible_navier_stokes_1d.jl
+++ b/src/equations/compressible_navier_stokes_1d.jl
@@ -148,7 +148,7 @@ end
 function flux(u, gradients, orientation::Integer,
               equations::CompressibleNavierStokesDiffusion1D)
     # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`.
-    rho, v1, _ = convert_transformed_to_primitive(u, equations)
+    _, v1, _ = convert_transformed_to_primitive(u, equations)
     # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, T)
     # either computed directly or reverse engineered from the gradient of the entropy variables
     # by way of the `convert_gradient_variables` function.
@@ -257,6 +257,12 @@ end
     return T
 end
 
+@inline function velocity(u, equations::CompressibleNavierStokesDiffusion1D)
+    rho = u[1]
+    v1 = u[2] / rho
+    return v1
+end
+
 @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip,
                                                                         <:Adiabatic})(flux_inner,
                                                                                       u_inner,
diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl
index 5708abb4e00..977dae18a43 100644
--- a/src/equations/compressible_navier_stokes_2d.jl
+++ b/src/equations/compressible_navier_stokes_2d.jl
@@ -148,7 +148,7 @@ end
 function flux(u, gradients, orientation::Integer,
               equations::CompressibleNavierStokesDiffusion2D)
     # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`.
-    rho, v1, v2, _ = convert_transformed_to_primitive(u, equations)
+    _, v1, v2, _ = convert_transformed_to_primitive(u, equations)
     # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, T)
     # either computed directly or reverse engineered from the gradient of the entropy variables
     # by way of the `convert_gradient_variables` function.
@@ -281,6 +281,13 @@ end
     return T
 end
 
+@inline function velocity(u, equations::CompressibleNavierStokesDiffusion2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    return SVector(v1, v2)
+end
+
 @inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion2D)
     # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v
 
diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl
index 9622882fc1a..551462e7595 100644
--- a/src/equations/compressible_navier_stokes_3d.jl
+++ b/src/equations/compressible_navier_stokes_3d.jl
@@ -148,7 +148,7 @@ end
 function flux(u, gradients, orientation::Integer,
               equations::CompressibleNavierStokesDiffusion3D)
     # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`.
-    rho, v1, v2, v3, _ = convert_transformed_to_primitive(u, equations)
+    _, v1, v2, v3, _ = convert_transformed_to_primitive(u, equations)
     # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, v3, T)
     # either computed directly or reverse engineered from the gradient of the entropy variables
     # by way of the `convert_gradient_variables` function.
@@ -309,6 +309,14 @@ end
     return T
 end
 
+@inline function velocity(u, equations::CompressibleNavierStokesDiffusion3D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    return SVector(v1, v2, v3)
+end
+
 @inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion3D)
     # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v
 
diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl
index 7de57a45178..c640be3c8b8 100644
--- a/test/test_parabolic_1d.jl
+++ b/test/test_parabolic_1d.jl
@@ -228,6 +228,29 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "TreeMesh1D: elixir_navierstokes_viscous_shock.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem",
+                                 "elixir_navierstokes_viscous_shock.jl"),
+                        l2=[
+                            0.00025762354103445303,
+                            0.0001433692781569829,
+                            0.00017369861968287976
+                        ],
+                        linf=[
+                            0.0016731940030498826,
+                            0.0010638575921477766,
+                            0.0011495207677434394
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 
 # Clean up afterwards: delete Trixi output directory
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index dc3e776154f..62142e1fc3b 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -760,6 +760,31 @@ end
     @test isapprox(drag_f, 1.5427441885921553, atol = 1e-13)
     @test isapprox(lift_f, 0.005621910087395724, atol = 1e-13)
 end
+
+@trixi_testset "P4estMesh2D: elixir_navierstokes_viscous_shock.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem",
+                                 "elixir_navierstokes_viscous_shock.jl"),
+                        l2=[
+                            0.0002576236264053728,
+                            0.00014336949098706463,
+                            7.189100338239794e-17,
+                            0.00017369905124642074
+                        ],
+                        linf=[
+                            0.0016731940983241156,
+                            0.0010638640749656147,
+                            5.59044079947959e-16,
+                            0.001149532023891009
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 
 # Clean up afterwards: delete Trixi.jl output directory
diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl
index 99ad05c7f70..9b7cb8a75ce 100644
--- a/test/test_parabolic_3d.jl
+++ b/test/test_parabolic_3d.jl
@@ -503,6 +503,33 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "P4estMesh3D: elixir_navierstokes_viscous_shock.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "p4est_3d_dgsem",
+                                 "elixir_navierstokes_viscous_shock.jl"),
+                        l2=[
+                            0.0002576235461250765,
+                            0.0001433693418567713,
+                            1.5583069105517042e-16,
+                            1.257551423107977e-16,
+                            0.00017369872990116004
+                        ],
+                        linf=[
+                            0.0016731930282756213,
+                            0.0010638586882356638,
+                            2.738015991633e-15,
+                            3.281831854493919e-15,
+                            0.0011495231318404686
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 
 # Clean up afterwards: delete Trixi.jl output directory