diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl
index a71750ff98..542a5b7039 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 e964d9c2b2..37a4f36210 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 4f4d4553a8..56a24238fc 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 da434579f7..c25a463118 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 2424ad0301..2278fe58eb 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 7875954c5f..0304ea16cf 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 3253e4410f..685ad6f8f3 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 7a6d19facd..a2c748ae6f 100644
--- a/src/equations/ideal_glm_mhd_2d.jl
+++ b/src/equations/ideal_glm_mhd_2d.jl
@@ -1119,6 +1119,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 e922a2e6fd..1fbf67134c 100644
--- a/src/equations/ideal_glm_mhd_3d.jl
+++ b/src/equations/ideal_glm_mhd_3d.jl
@@ -1010,6 +1010,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 571d4723f6..98f0f7f1dc 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 4ecaf3b6e1..8b19d19a21 100644
--- a/src/equations/shallow_water_2d.jl
+++ b/src/equations/shallow_water_2d.jl
@@ -283,7 +283,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)
 """
@@ -840,6 +840,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
@@ -903,11 +909,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 d507cded95..58e1fa6ef2 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)) ==
@@ -1165,6 +1177,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_ll,
@@ -1183,6 +1196,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
@@ -1221,6 +1235,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
@@ -1283,6 +1298,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_ll,
@@ -1301,6 +1317,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
@@ -1319,6 +1336,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
@@ -2086,6 +2104,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)) ==
@@ -2101,6 +2120,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)) ==
@@ -2113,6 +2133,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
@@ -2240,6 +2261,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_ll,
@@ -2285,6 +2307,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 5831122ffe..bccdcf8faa 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