diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl
index a71750ff98c..542a5b7039f 100644
--- a/src/equations/compressible_euler_1d.jl
+++ b/src/equations/compressible_euler_1d.jl
@@ -412,10 +412,10 @@ See also
     return SVector(f1, f2, f3)
 end
 
-# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that 
-# the normal component is incorporated into the numerical flux. 
-# 
-# See `flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})` for a 
+# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that
+# the normal component is incorporated into the numerical flux.
+#
+# See `flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})` for a
 # similar implementation.
 @inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector,
                               equations::CompressibleEulerEquations1D)
@@ -943,6 +943,12 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::CompressibleEulerEquations1D)
+    rho = u[1]
+    v1 = u[2] / rho
+    return v1
+end
+
 @inline function pressure(u, equations::CompressibleEulerEquations1D)
     rho, rho_v1, rho_e = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2) / rho)
diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
index e964d9c2b21..37a4f362101 100644
--- a/src/equations/compressible_euler_2d.jl
+++ b/src/equations/compressible_euler_2d.jl
@@ -1957,6 +1957,19 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::CompressibleEulerEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    return SVector(v1, v2)
+end
+
+@inline function velocity(u, orientation::Int, equations::CompressibleEulerEquations2D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
 @inline function pressure(u, equations::CompressibleEulerEquations2D)
     rho, rho_v1, rho_v2, rho_e = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2) / rho)
diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl
index 4f4d4553a8f..56a24238fcc 100644
--- a/src/equations/compressible_euler_3d.jl
+++ b/src/equations/compressible_euler_3d.jl
@@ -1692,6 +1692,20 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::CompressibleEulerEquations3D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    return SVector(v1, v2, v3)
+end
+
+@inline function velocity(u, orientation::Int, equations::CompressibleEulerEquations3D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
 @inline function pressure(u, equations::CompressibleEulerEquations3D)
     rho, rho_v1, rho_v2, rho_v3, rho_e = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho)
diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl
index da434579f76..c25a4631185 100644
--- a/src/equations/compressible_euler_multicomponent_1d.jl
+++ b/src/equations/compressible_euler_multicomponent_1d.jl
@@ -618,4 +618,10 @@ end
     return SVector{ncomponents(equations), real(equations)}(u[i + 2] * v
                                                             for i in eachcomponent(equations))
 end
+
+@inline function velocity(u, equations::CompressibleEulerMulticomponentEquations1D)
+    rho = density(u, equations)
+    v1 = u[1] / rho
+    return v1
+end
 end # @muladd
diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl
index 2424ad0301c..2278fe58eb3 100644
--- a/src/equations/compressible_euler_multicomponent_2d.jl
+++ b/src/equations/compressible_euler_multicomponent_2d.jl
@@ -825,4 +825,18 @@ end
     return SVector{ncomponents(equations), real(equations)}(u[i + 3] * v
                                                             for i in eachcomponent(equations))
 end
+
+@inline function velocity(u, equations::CompressibleEulerMulticomponentEquations2D)
+    rho = density(u, equations)
+    v1 = u[1] / rho
+    v2 = u[2] / rho
+    return SVector(v1, v2)
+end
+
+@inline function velocity(u, orientation::Int,
+                          equations::CompressibleEulerMulticomponentEquations2D)
+    rho = density(u, equations)
+    v = u[orientation] / rho
+    return v
+end
 end # @muladd
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index 7875954c5f0..0304ea16cfe 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -307,6 +307,37 @@ The inverse conversion is performed by [`cons2prim`](@ref).
 """
 function prim2cons end
 
+"""
+    velocity(u, equations)
+
+Return the velocity vector corresponding to the equations, e.g., fluid velocity for
+Euler's equations. The velocity in certain orientation or normal direction (scalar) can be computed
+with `velocity(u, orientation, equations)` or `velocity(u, normal_direction, equations)`
+respectively. The `velocity(u, normal_direction, equations)` function calls
+`velocity(u, equations)` to compute the velocity vector and then normal vector, thus allowing
+a general function to be written for the AbstractEquations type. However, the
+`velocity(u, orientation, equations)` is written for each equation separately to ensure
+only the velocity in the desired direction (orientation) is computed.
+`u` is a vector of the conserved variables at a single node, i.e., a vector
+of the correct length `nvariables(equations)`.
+"""
+function velocity end
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::AbstractEquations{2})
+    vel = velocity(u, equations)
+    v = vel[1] * normal_direction[1] + vel[2] * normal_direction[2]
+    return v
+end
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::AbstractEquations{3})
+    vel = velocity(u, equations)
+    v = vel[1] * normal_direction[1] + vel[2] * normal_direction[2] +
+        vel[3] * normal_direction[3]
+    return v
+end
+
 """
     entropy(u, equations)
 
diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl
index 3253e4410f8..685ad6f8f3f 100644
--- a/src/equations/ideal_glm_mhd_1d.jl
+++ b/src/equations/ideal_glm_mhd_1d.jl
@@ -326,7 +326,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
         Sdiff = SsR - SsL
 
         # Compute HLL values for vStar, BStar
-        # These correspond to eq. (28) and (30) from the referenced paper 
+        # These correspond to eq. (28) and (30) from the referenced paper
         # and the classic HLL intermediate state given by (2)
         rho_HLL = (SsR * rho_rr - SsL * rho_ll - (f_rr[1] - f_ll[1])) / Sdiff
 
@@ -394,7 +394,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
                          B1_rr * (B1_rr * v1_rr + B2_rr * v2_rr + B3_rr * v3_rr))) /
                        SdiffStar # (23)
 
-            # Classic HLLC update (32)           
+            # Classic HLLC update (32)
             f1 = f_rr[1] + SsR * (densStar - u_rr[1])
             f2 = f_rr[2] + SsR * (mom_1_Star - u_rr[2])
             f3 = f_rr[3] + SsR * (mom_2_Star - u_rr[3])
@@ -555,6 +555,20 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::IdealGlmMhdEquations1D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    return SVector(v1, v2, v3)
+end
+
+@inline function velocity(u, orientation::Int, equations::IdealGlmMhdEquations1D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
 @inline function pressure(u, equations::IdealGlmMhdEquations1D)
     rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
index a43c6b0e619..9383fdbe961 100644
--- a/src/equations/ideal_glm_mhd_2d.jl
+++ b/src/equations/ideal_glm_mhd_2d.jl
@@ -1111,6 +1111,20 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::IdealGlmMhdEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    return SVector(v1, v2, v3)
+end
+
+@inline function velocity(u, orientation::Int, equations::IdealGlmMhdEquations2D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
 @inline function pressure(u, equations::IdealGlmMhdEquations2D)
     rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl
index 9bb05cb3f89..8d7d0461398 100644
--- a/src/equations/ideal_glm_mhd_3d.jl
+++ b/src/equations/ideal_glm_mhd_3d.jl
@@ -1001,6 +1001,20 @@ end
     return u[1]
 end
 
+@inline function velocity(u, equations::IdealGlmMhdEquations3D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    return SVector(v1, v2, v3)
+end
+
+@inline function velocity(u, orientation::Int, equations::IdealGlmMhdEquations3D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
 @inline function pressure(u, equations::IdealGlmMhdEquations3D)
     rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl
index 571d4723f6f..98f0f7f1dc1 100644
--- a/src/equations/polytropic_euler_2d.jl
+++ b/src/equations/polytropic_euler_2d.jl
@@ -389,6 +389,19 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::PolytropicEulerEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    return SVector(v1, v2)
+end
+
+@inline function velocity(u, orientation::Int, equations::PolytropicEulerEquations2D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
 @inline function pressure(u, equations::PolytropicEulerEquations2D)
     rho, rho_v1, rho_v2 = u
     p = equations.kappa * rho^equations.gamma
diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl
index 30eca095aa4..7ba36d6ec92 100644
--- a/src/equations/shallow_water_2d.jl
+++ b/src/equations/shallow_water_2d.jl
@@ -282,7 +282,7 @@ Further details are available in the papers:
   shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry
   [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036)
 - Patrick Ersing, Andrew R. Winters (2023)
-  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
+  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on
   curvilinear meshes
   [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
 """
@@ -834,6 +834,12 @@ end
     return SVector(v1, v2)
 end
 
+@inline function velocity(u, orientation::Int, equations::ShallowWaterEquations2D)
+    h = u[1]
+    v = u[orientation + 1] / h
+    return v
+end
+
 # Convert conservative variables to primitive
 @inline function cons2prim(u, equations::ShallowWaterEquations2D)
     h, _, _, b = u
@@ -897,11 +903,11 @@ end
                        equations::ShallowWaterEquations2D)
 
 Calculate Roe-averaged velocity `v_roe` and wavespeed `c_roe = sqrt{g * h_roe}` depending on direction.
-See for instance equation (62) in 
+See for instance equation (62) in
 - Paul A. Ullrich, Christiane Jablonowski, and Bram van Leer (2010)
   High-order finite-volume methods for the shallow-water equations on the sphere
   [DOI: 10.1016/j.jcp.2010.04.044](https://doi.org/10.1016/j.jcp.2010.04.044)
-Or [this slides](https://faculty.washington.edu/rjl/classes/am574w2011/slides/am574lecture20nup3.pdf), 
+Or [this slides](https://faculty.washington.edu/rjl/classes/am574w2011/slides/am574lecture20nup3.pdf),
 slides 8 and 9.
 """
 @inline function calc_wavespeed_roe(u_ll, u_rr, orientation::Integer,
diff --git a/test/test_type.jl b/test/test_type.jl
index 62569efe0ec..92274bf4248 100644
--- a/test/test_type.jl
+++ b/test/test_type.jl
@@ -167,6 +167,7 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test eltype(@inferred cons2entropy(u, equations)) == RealT
             @test eltype(@inferred entropy2cons(u, equations)) == RealT
             @test typeof(@inferred density(u, equations)) == RealT
+            @test typeof(@inferred velocity(u, equations)) == RealT
             @test typeof(@inferred pressure(u, equations)) == RealT
             @test typeof(@inferred density_pressure(u, equations)) == RealT
             @test typeof(@inferred entropy(cons, equations)) == RealT
@@ -218,6 +219,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                 end
             end
 
+            @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_shima_etal(u_ll, u_rr, normal_direction, equations)) ==
                   RealT
@@ -257,6 +259,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                   RealT
 
             for orientation in orientations
+                @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux_shima_etal(u_ll, u_rr, orientation, equations)) ==
                       RealT
@@ -300,6 +303,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                       RealT
             end
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
             @test eltype(@inferred cons2prim(u, equations)) == RealT
             @test eltype(@inferred prim2cons(u, equations)) == RealT
@@ -367,6 +371,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                 end
             end
 
+            @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_shima_etal(u_ll, u_rr, normal_direction, equations)) ==
                   RealT
@@ -392,6 +397,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                           equations)) == RealT
 
             for orientation in orientations
+                @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux_shima_etal(u_ll, u_rr, orientation, equations)) ==
                       RealT
@@ -421,6 +427,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                               equations)) == RealT
             end
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
             @test eltype(@inferred cons2prim(u, equations)) == RealT
             @test eltype(@inferred prim2cons(u, equations)) == RealT
@@ -456,6 +463,7 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) ==
                   RealT
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred flux(u, orientation, equations)) == RealT
             @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation,
                                                       equations)) == RealT
@@ -499,11 +507,13 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) ==
                   RealT
 
+            @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction,
                                                 equations)) == RealT
 
             for orientation in orientations
+                @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation,
                                                           equations)) == RealT
@@ -515,6 +525,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                       RealT
             end
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
             @test eltype(@inferred cons2prim(u, equations)) == RealT
             @test eltype(@inferred prim2cons(u, equations)) == RealT
@@ -1094,6 +1105,7 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) ==
                   RealT
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred flux(u, orientation, equations)) == RealT
             @test eltype(@inferred flux_hllc(u_ll, u_rr, orientation, equations)) == RealT
             @test eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, equations)) ==
@@ -1164,6 +1176,7 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) ==
                   RealT
 
+            @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr,
                                                                normal_direction,
@@ -1181,6 +1194,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                           equations)) == RealT
 
             for orientation in orientations
+                @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, orientation,
                                                                    equations)) == RealT
@@ -1219,6 +1233,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                               equations)) == RealT
             end
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
             @test eltype(@inferred cons2prim(u, equations)) == RealT
             @test eltype(@inferred prim2cons(u, equations)) == RealT
@@ -1279,6 +1294,7 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) ==
                   RealT
 
+            @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr,
                                                                normal_direction,
@@ -1296,6 +1312,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                           equations)) == RealT
 
             for orientation in orientations
+                @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, orientation,
                                                                    equations)) == RealT
@@ -1314,6 +1331,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                               equations)) == RealT
             end
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
             @test eltype(@inferred cons2prim(u, equations)) == RealT
             @test eltype(@inferred prim2cons(u, equations)) == RealT
@@ -2081,6 +2099,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                 @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) ==
                       RealT
 
+                @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
                 @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
                 @test eltype(@inferred flux_winters_etal(u_ll, u_rr, normal_direction,
                                                          equations)) ==
@@ -2096,6 +2115,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                       RealT
 
                 for orientation in orientations
+                    @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                     @test eltype(@inferred flux(u, orientation, equations)) == RealT
                     @test eltype(@inferred flux_winters_etal(u_ll, u_rr, orientation,
                                                              equations)) ==
@@ -2108,6 +2128,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                           RealT
                 end
 
+                @test eltype(@inferred velocity(u, equations)) == RealT
                 @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
                 @test eltype(@inferred cons2prim(u, equations)) == RealT
                 @test eltype(@inferred prim2cons(u, equations)) == RealT
@@ -2234,6 +2255,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                                 surface_flux_function,
                                                                 equations)) == RealT
 
+            @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
                                                                          normal_direction,
@@ -2276,6 +2298,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                           RealT
                 end
 
+                @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
                                                                              orientation,
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 5831122ffe2..bccdcf8faaa 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1737,6 +1737,123 @@ end
     @test isapprox(mu_control(u, equations_parabolic, T_ref, R_specific, C, mu_ref),
                    1.803e-5, atol = 5e-8)
 end
+
+# Velocity functions are present in many equations and are tested here
+@testset "Velocity functions for different equations" begin
+    gamma = 1.4
+    rho = pi * pi
+    pres = sqrt(pi)
+    v1, v2, v3 = pi, exp(1.0), exp(pi) # use pi, exp to test with non-trivial numbers
+    v_vector = SVector(v1, v2, v3)
+    normal_direction_2d = SVector(pi^2, pi^3)
+    normal_direction_3d = SVector(normal_direction_2d..., pi^4)
+    v_normal_1d = v1 * normal_direction_2d[1]
+    v_normal_2d = v1 * normal_direction_2d[1] + v2 * normal_direction_2d[2]
+    v_normal_3d = v_normal_2d + v3 * normal_direction_3d[3]
+
+    equations_euler_1d = CompressibleEulerEquations1D(gamma)
+    u = prim2cons(SVector(rho, v1, pres), equations_euler_1d)
+    @test isapprox(velocity(u, equations_euler_1d), v1)
+
+    equations_euler_2d = CompressibleEulerEquations2D(gamma)
+    u = prim2cons(SVector(rho, v1, v2, pres), equations_euler_2d)
+    @test isapprox(velocity(u, equations_euler_2d), SVector(v1, v2))
+    @test isapprox(velocity(u, normal_direction_2d, equations_euler_2d), v_normal_2d)
+    for orientation in 1:2
+        @test isapprox(velocity(u, orientation, equations_euler_2d),
+                       v_vector[orientation])
+    end
+
+    equations_euler_3d = CompressibleEulerEquations3D(gamma)
+    u = prim2cons(SVector(rho, v1, v2, v3, pres), equations_euler_3d)
+    @test isapprox(velocity(u, equations_euler_3d), SVector(v1, v2, v3))
+    @test isapprox(velocity(u, normal_direction_3d, equations_euler_3d), v_normal_3d)
+    for orientation in 1:3
+        @test isapprox(velocity(u, orientation, equations_euler_3d),
+                       v_vector[orientation])
+    end
+
+    rho1, rho2 = rho, rho * pi # use pi to test with non-trivial numbers
+    gammas = (gamma, exp(gamma))
+    gas_constants = (0.387, 1.678) # Standard numbers + 0.1
+
+    equations_multi_euler_1d = CompressibleEulerMulticomponentEquations1D(; gammas,
+                                                                          gas_constants)
+    u = prim2cons(SVector(v1, pres, rho1, rho2), equations_multi_euler_1d)
+    @test isapprox(velocity(u, equations_multi_euler_1d), v1)
+
+    equations_multi_euler_2d = CompressibleEulerMulticomponentEquations2D(; gammas,
+                                                                          gas_constants)
+    u = prim2cons(SVector(v1, v2, pres, rho1, rho2), equations_multi_euler_2d)
+    @test isapprox(velocity(u, equations_multi_euler_2d), SVector(v1, v2))
+    @test isapprox(velocity(u, normal_direction_2d, equations_multi_euler_2d),
+                   v_normal_2d)
+    for orientation in 1:2
+        @test isapprox(velocity(u, orientation, equations_multi_euler_2d),
+                       v_vector[orientation])
+    end
+
+    kappa = 0.1 * pi # pi for non-trivial test
+    equations_polytropic = PolytropicEulerEquations2D(gamma, kappa)
+    u = prim2cons(SVector(rho, v1, v2), equations_polytropic)
+    @test isapprox(velocity(u, equations_polytropic), SVector(v1, v2))
+    equations_polytropic = CompressibleEulerMulticomponentEquations2D(; gammas,
+                                                                      gas_constants)
+    u = prim2cons(SVector(v1, v2, pres, rho1, rho2), equations_polytropic)
+    @test isapprox(velocity(u, equations_polytropic), SVector(v1, v2))
+    @test isapprox(velocity(u, normal_direction_2d, equations_polytropic), v_normal_2d)
+    for orientation in 1:2
+        @test isapprox(velocity(u, orientation, equations_polytropic),
+                       v_vector[orientation])
+    end
+
+    B1, B2, B3 = pi^3, pi^4, pi^5
+    equations_ideal_mhd_1d = IdealGlmMhdEquations1D(gamma)
+    u = prim2cons(SVector(rho, v1, v2, v3, pres, B1, B2, B3), equations_ideal_mhd_1d)
+    @test isapprox(velocity(u, equations_ideal_mhd_1d), SVector(v1, v2, v3))
+    for orientation in 1:3
+        @test isapprox(velocity(u, orientation, equations_ideal_mhd_1d),
+                       v_vector[orientation])
+    end
+
+    psi = exp(0.1)
+    equations_ideal_mhd_2d = IdealGlmMhdEquations2D(gamma)
+    u = prim2cons(SVector(rho, v1, v2, v3, pres, B1, B2, B3, psi),
+                  equations_ideal_mhd_2d)
+    @test isapprox(velocity(u, equations_ideal_mhd_2d), SVector(v1, v2, v3))
+    @test isapprox(velocity(u, normal_direction_2d, equations_ideal_mhd_2d),
+                   v_normal_2d)
+    for orientation in 1:3
+        @test isapprox(velocity(u, orientation, equations_ideal_mhd_2d),
+                       v_vector[orientation])
+    end
+
+    equations_ideal_mhd_3d = IdealGlmMhdEquations3D(gamma)
+    u = prim2cons(SVector(rho, v1, v2, v3, pres, B1, B2, B3, psi),
+                  equations_ideal_mhd_3d)
+    @test isapprox(velocity(u, equations_ideal_mhd_3d), SVector(v1, v2, v3))
+    @test isapprox(velocity(u, normal_direction_3d, equations_ideal_mhd_3d),
+                   v_normal_3d)
+    for orientation in 1:3
+        @test isapprox(velocity(u, orientation, equations_ideal_mhd_3d),
+                       v_vector[orientation])
+    end
+
+    H, b = exp(pi), exp(pi^2)
+    gravity_constant, H0 = 9.91, 0.1 # Standard numbers + 0.1
+    shallow_water_1d = ShallowWaterEquations1D(; gravity_constant, H0)
+    u = prim2cons(SVector(H, v1, b), shallow_water_1d)
+    @test isapprox(velocity(u, shallow_water_1d), v1)
+
+    shallow_water_2d = ShallowWaterEquations2D(; gravity_constant, H0)
+    u = prim2cons(SVector(H, v1, v2, b), shallow_water_2d)
+    @test isapprox(velocity(u, shallow_water_2d), SVector(v1, v2))
+    @test isapprox(velocity(u, normal_direction_2d, shallow_water_2d), v_normal_2d)
+    for orientation in 1:2
+        @test isapprox(velocity(u, orientation, shallow_water_2d),
+                       v_vector[orientation])
+    end
+end
 end
 
 end #module