From 217788acc52c9cf333be15b0f5798586809fa407 Mon Sep 17 00:00:00 2001
From: Arpit-Babbar <arpitbabbar@gmail.com>
Date: Tue, 17 Sep 2024 06:42:55 +0200
Subject: [PATCH 1/9] Add velocity functions for different equations

---
 src/equations/compressible_euler_1d.jl        | 14 ++++++++---
 src/equations/compressible_euler_2d.jl        | 22 ++++++++++++++++
 src/equations/compressible_euler_3d.jl        | 24 ++++++++++++++++++
 .../compressible_euler_multicomponent_1d.jl   |  8 ++++++
 .../compressible_euler_multicomponent_2d.jl   | 23 +++++++++++++++++
 src/equations/ideal_glm_mhd_1d.jl             | 10 ++++++--
 src/equations/ideal_glm_mhd_2d.jl             | 22 ++++++++++++++++
 src/equations/ideal_glm_mhd_3d.jl             | 25 +++++++++++++++++++
 src/equations/polytropic_euler_2d.jl          | 13 ++++++++++
 src/equations/shallow_water_2d.jl             | 21 +++++++++++++---
 test/test_type.jl                             | 24 ++++++++++++++++++
 11 files changed, 197 insertions(+), 9 deletions(-)

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..7884e65d6b 100644
--- a/src/equations/compressible_euler_2d.jl
+++ b/src/equations/compressible_euler_2d.jl
@@ -1957,6 +1957,28 @@ 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 velocity(u, normal_direction::AbstractVector,
+                          equations::CompressibleEulerEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2]
+    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..79ade4e7ae 100644
--- a/src/equations/compressible_euler_3d.jl
+++ b/src/equations/compressible_euler_3d.jl
@@ -1692,6 +1692,30 @@ 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 velocity(u, normal_direction::AbstractVector,
+                          equations::CompressibleEulerEquations3D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2] + v3 * normal_direction[3]
+    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..ec13b1f84e 100644
--- a/src/equations/compressible_euler_multicomponent_1d.jl
+++ b/src/equations/compressible_euler_multicomponent_1d.jl
@@ -618,4 +618,12 @@ 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..6e2099b2ab 100644
--- a/src/equations/compressible_euler_multicomponent_2d.jl
+++ b/src/equations/compressible_euler_multicomponent_2d.jl
@@ -825,4 +825,27 @@ 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
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::CompressibleEulerMulticomponentEquations2D)
+                          rho = density(u, equations)
+    v1 = u[1] / rho
+    v2 = u[2] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2]
+    return v
+end
+
 end # @muladd
diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl
index 3253e4410f..cbfe37daa7 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,12 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::IdealGlmMhdEquations1D)
+    rho = u[1]
+    v1 = u[2] / rho
+    return v1
+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..801f89afab 100644
--- a/src/equations/ideal_glm_mhd_2d.jl
+++ b/src/equations/ideal_glm_mhd_2d.jl
@@ -1119,6 +1119,28 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::IdealGlmMhdEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    return SVector(v1, v2)
+end
+
+@inline function velocity(u, orientation::Int, equations::IdealGlmMhdEquations2D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::IdealGlmMhdEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2]
+    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..afbddbb78b 100644
--- a/src/equations/ideal_glm_mhd_3d.jl
+++ b/src/equations/ideal_glm_mhd_3d.jl
@@ -1010,6 +1010,31 @@ 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 velocity(u, normal_direction::AbstractVector,
+                          equations::IdealGlmMhdEquations3D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2] + v3 * normal_direction[3]
+    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..e060b389ff 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, 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..3f5cc082dd 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,21 @@ 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
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::ShallowWaterEquations2D)
+    h = u[1]
+    v1 = u[2] / h
+    v2 = u[3] / h
+    v = v1 * normal_direction[1] + v2 * normal_direction[2]
+    return v
+end
+
 # Convert conservative variables to primitive
 @inline function cons2prim(u, equations::ShallowWaterEquations2D)
     h, _, _, b = u
@@ -903,11 +918,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..9787bdd956 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
@@ -1425,6 +1443,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_nonconservative_powell(u_ll, u_rr, orientation,
                                                                    equations)) ==
@@ -2086,6 +2105,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 +2121,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 +2134,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 +2262,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 +2308,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,

From 9d2a0eb3593851fac647515b92de5fc976b5060a Mon Sep 17 00:00:00 2001
From: Arpit Babbar <arpitbabbar@gmail.com>
Date: Mon, 23 Sep 2024 15:36:42 +0200
Subject: [PATCH 2/9] Formatter

---
 src/equations/compressible_euler_multicomponent_1d.jl | 2 --
 src/equations/compressible_euler_multicomponent_2d.jl | 6 +++---
 src/equations/ideal_glm_mhd_3d.jl                     | 1 -
 3 files changed, 3 insertions(+), 6 deletions(-)

diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl
index ec13b1f84e..c25a463118 100644
--- a/src/equations/compressible_euler_multicomponent_1d.jl
+++ b/src/equations/compressible_euler_multicomponent_1d.jl
@@ -624,6 +624,4 @@ end
     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 6e2099b2ab..65098f307b 100644
--- a/src/equations/compressible_euler_multicomponent_2d.jl
+++ b/src/equations/compressible_euler_multicomponent_2d.jl
@@ -833,7 +833,8 @@ end
     return SVector(v1, v2)
 end
 
-@inline function velocity(u, orientation::Int, equations::CompressibleEulerMulticomponentEquations2D)
+@inline function velocity(u, orientation::Int,
+                          equations::CompressibleEulerMulticomponentEquations2D)
     rho = density(u, equations)
     v = u[orientation] / rho
     return v
@@ -841,11 +842,10 @@ end
 
 @inline function velocity(u, normal_direction::AbstractVector,
                           equations::CompressibleEulerMulticomponentEquations2D)
-                          rho = density(u, equations)
+    rho = density(u, equations)
     v1 = u[1] / rho
     v2 = u[2] / rho
     v = v1 * normal_direction[1] + v2 * normal_direction[2]
     return v
 end
-
 end # @muladd
diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl
index afbddbb78b..5bffd23c84 100644
--- a/src/equations/ideal_glm_mhd_3d.jl
+++ b/src/equations/ideal_glm_mhd_3d.jl
@@ -1034,7 +1034,6 @@ end
     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

From 7386197b34cabf4e739b3acc4c78c5116472ad34 Mon Sep 17 00:00:00 2001
From: Arpit Babbar <arpitbabbar@gmail.com>
Date: Mon, 23 Sep 2024 16:15:59 +0200
Subject: [PATCH 3/9] Fix CI

---
 src/equations/polytropic_euler_2d.jl | 10 +++++++++-
 test/test_type.jl                    |  1 -
 2 files changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl
index e060b389ff..fca89d1195 100644
--- a/src/equations/polytropic_euler_2d.jl
+++ b/src/equations/polytropic_euler_2d.jl
@@ -396,12 +396,20 @@ end
     return SVector(v1, v2)
 end
 
-@inline function velocity(u, orientation, equations::PolytropicEulerEquations2D)
+@inline function velocity(u, orientation::Int, equations::PolytropicEulerEquations2D)
     rho = u[1]
     v = u[orientation + 1] / rho
     return v
 end
 
+@inline function velocity(u, normal_direction::AbstractVector, equations::PolytropicEulerEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2]
+    return v
+end
+
 @inline function pressure(u, equations::PolytropicEulerEquations2D)
     rho, rho_v1, rho_v2 = u
     p = equations.kappa * rho^equations.gamma
diff --git a/test/test_type.jl b/test/test_type.jl
index 9787bdd956..58e1fa6ef2 100644
--- a/test/test_type.jl
+++ b/test/test_type.jl
@@ -1443,7 +1443,6 @@ 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_nonconservative_powell(u_ll, u_rr, orientation,
                                                                    equations)) ==

From 35b776dceca1124a40e2647874b343b3c809228d Mon Sep 17 00:00:00 2001
From: Arpit Babbar <arpitbabbar@gmail.com>
Date: Mon, 23 Sep 2024 16:17:53 +0200
Subject: [PATCH 4/9] Format

---
 src/equations/polytropic_euler_2d.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl
index fca89d1195..baa05de820 100644
--- a/src/equations/polytropic_euler_2d.jl
+++ b/src/equations/polytropic_euler_2d.jl
@@ -402,7 +402,8 @@ end
     return v
 end
 
-@inline function velocity(u, normal_direction::AbstractVector, equations::PolytropicEulerEquations2D)
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::PolytropicEulerEquations2D)
     rho = u[1]
     v1 = u[2] / rho
     v2 = u[3] / rho

From 893accd6f0017b614bdf4e1590d247cfdd01937f Mon Sep 17 00:00:00 2001
From: Arpit Babbar <arpitbabbar@gmail.com>
Date: Tue, 24 Sep 2024 16:56:27 +0200
Subject: [PATCH 5/9] Export velocity, add docs

---
 src/Trixi.jl               |  2 +-
 src/equations/equations.jl | 13 +++++++++++++
 2 files changed, 14 insertions(+), 1 deletion(-)

diff --git a/src/Trixi.jl b/src/Trixi.jl
index d9b9245918..ec32c52ed9 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -219,7 +219,7 @@ export initial_condition_eoc_test_coupled_euler_gravity,
 
 export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean,
        cons2entropy, entropy2cons
-export density, pressure, density_pressure, velocity, global_mean_vars,
+export density, velocity, pressure, density_pressure, velocity, global_mean_vars,
        equilibrium_distribution, waterheight_pressure
 export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic,
        cross_helicity,
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index 7875954c5f..a95e619a3b 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -307,6 +307,19 @@ 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.
+
+`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
+
 """
     entropy(u, equations)
 

From 6f2ae514a7a8737ca1ddc498f589572339dcc9a2 Mon Sep 17 00:00:00 2001
From: Arpit Babbar <arpitbabbar@gmail.com>
Date: Wed, 25 Sep 2024 16:58:45 +0200
Subject: [PATCH 6/9] Export velocity, add docstring, add unit tests, make MHD
 velocity like 3D for all dims

---
 src/equations/ideal_glm_mhd_1d.jl |  20 ++++-
 src/equations/ideal_glm_mhd_2d.jl |   6 +-
 test/test_unit.jl                 | 120 +++++++++++++++++++++++++++++-
 3 files changed, 142 insertions(+), 4 deletions(-)

diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl
index cbfe37daa7..5160355804 100644
--- a/src/equations/ideal_glm_mhd_1d.jl
+++ b/src/equations/ideal_glm_mhd_1d.jl
@@ -558,7 +558,25 @@ end
 @inline function velocity(u, equations::IdealGlmMhdEquations1D)
     rho = u[1]
     v1 = u[2] / rho
-    return v1
+    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 velocity(u, normal_direction::AbstractVector,
+                          equations::IdealGlmMhdEquations1D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2] + v3 * normal_direction[3]
+    return v
 end
 
 @inline function pressure(u, equations::IdealGlmMhdEquations1D)
diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
index 801f89afab..614e6a4cd7 100644
--- a/src/equations/ideal_glm_mhd_2d.jl
+++ b/src/equations/ideal_glm_mhd_2d.jl
@@ -1123,7 +1123,8 @@ end
     rho = u[1]
     v1 = u[2] / rho
     v2 = u[3] / rho
-    return SVector(v1, v2)
+    v3 = u[4] / rho
+    return SVector(v1, v2, v3)
 end
 
 @inline function velocity(u, orientation::Int, equations::IdealGlmMhdEquations2D)
@@ -1137,7 +1138,8 @@ end
     rho = u[1]
     v1 = u[2] / rho
     v2 = u[3] / rho
-    v = v1 * normal_direction[1] + v2 * normal_direction[2]
+    v3 = u[4] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2] + v3 * normal_direction[3]
     return v
 end
 
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 70e2e2ed10..66efdb0661 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1295,7 +1295,7 @@ end
 
     # Define wrapper function for pressure in order to call default implementation
     function pressure_test(u, equations)
-        return pressure(u, equations)
+        return pres(u, equations)
     end
 
     @test Trixi.gradient_conservative(pressure_test, u, equations) ≈
@@ -1737,6 +1737,124 @@ 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.4f0
+    rho = pi * pi
+    pres = sqrt(pi)
+    v1, v2, v3 = pi, exp(1.0f0), 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_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))
+    @test isapprox(velocity(u, normal_direction_3d, equations_ideal_mhd_1d),
+                   v_normal_3d)
+    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_3d, equations_ideal_mhd_2d),
+                   v_normal_3d)
+    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

From 7428094dff07e4066dbad85d53fa133f7663146c Mon Sep 17 00:00:00 2001
From: Arpit Babbar <arpitbabbar@gmail.com>
Date: Wed, 25 Sep 2024 17:00:10 +0200
Subject: [PATCH 7/9] Refractor bug fixgit push

---
 test/test_unit.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_unit.jl b/test/test_unit.jl
index 66efdb0661..9daee40a06 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1295,7 +1295,7 @@ end
 
     # Define wrapper function for pressure in order to call default implementation
     function pressure_test(u, equations)
-        return pres(u, equations)
+        return pressure(u, equations)
     end
 
     @test Trixi.gradient_conservative(pressure_test, u, equations) ≈

From cdf26532bef0c531f3fdbaeab169788442d5782c Mon Sep 17 00:00:00 2001
From: Arpit Babbar <arpitbabbar@gmail.com>
Date: Thu, 26 Sep 2024 07:43:36 +0200
Subject: [PATCH 8/9] Changes as per review

---
 src/Trixi.jl                      |  2 +-
 src/equations/ideal_glm_mhd_1d.jl |  4 +---
 src/equations/ideal_glm_mhd_2d.jl |  3 +--
 test/test_unit.jl                 | 13 +++++++------
 4 files changed, 10 insertions(+), 12 deletions(-)

diff --git a/src/Trixi.jl b/src/Trixi.jl
index ec32c52ed9..d9b9245918 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -219,7 +219,7 @@ export initial_condition_eoc_test_coupled_euler_gravity,
 
 export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean,
        cons2entropy, entropy2cons
-export density, velocity, pressure, density_pressure, velocity, global_mean_vars,
+export density, pressure, density_pressure, velocity, global_mean_vars,
        equilibrium_distribution, waterheight_pressure
 export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic,
        cross_helicity,
diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl
index 5160355804..24dd03a0d6 100644
--- a/src/equations/ideal_glm_mhd_1d.jl
+++ b/src/equations/ideal_glm_mhd_1d.jl
@@ -573,9 +573,7 @@ end
                           equations::IdealGlmMhdEquations1D)
     rho = u[1]
     v1 = u[2] / rho
-    v2 = u[3] / rho
-    v3 = u[4] / rho
-    v = v1 * normal_direction[1] + v2 * normal_direction[2] + v3 * normal_direction[3]
+    v = v1 * normal_direction[1]
     return v
 end
 
diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
index 614e6a4cd7..8d30920869 100644
--- a/src/equations/ideal_glm_mhd_2d.jl
+++ b/src/equations/ideal_glm_mhd_2d.jl
@@ -1138,8 +1138,7 @@ end
     rho = u[1]
     v1 = u[2] / rho
     v2 = u[3] / rho
-    v3 = u[4] / rho
-    v = v1 * normal_direction[1] + v2 * normal_direction[2] + v3 * normal_direction[3]
+    v = v1 * normal_direction[1] + v2 * normal_direction[2]
     return v
 end
 
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 9daee40a06..7bf2448247 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1740,13 +1740,14 @@ end
 
 # Velocity functions are present in many equations and are tested here
 @testset "Velocity functions for different equations" begin
-    gamma = 1.4f0
+    gamma = 1.4
     rho = pi * pi
     pres = sqrt(pi)
-    v1, v2, v3 = pi, exp(1.0f0), exp(pi) # use pi, exp to test with non-trivial numbers
+    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]
 
@@ -1810,8 +1811,8 @@ end
     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))
-    @test isapprox(velocity(u, normal_direction_3d, equations_ideal_mhd_1d),
-                   v_normal_3d)
+    @test isapprox(velocity(u, SVector(normal_direction_2d[1]), equations_ideal_mhd_1d),
+                   v_normal_1d)
     for orientation in 1:3
         @test isapprox(velocity(u, orientation, equations_ideal_mhd_1d),
                        v_vector[orientation])
@@ -1822,8 +1823,8 @@ end
     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_3d, equations_ideal_mhd_2d),
-                   v_normal_3d)
+    @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])

From d29df2b4fa56c6349530f8a26f75488f9dfe78dc Mon Sep 17 00:00:00 2001
From: Arpit Babbar <arpitbabbar@gmail.com>
Date: Sun, 29 Sep 2024 20:56:38 +0200
Subject: [PATCH 9/9] General velocity normal function

---
 src/equations/compressible_euler_2d.jl        |  9 --------
 src/equations/compressible_euler_3d.jl        | 10 ---------
 .../compressible_euler_multicomponent_2d.jl   |  9 --------
 src/equations/equations.jl                    | 22 +++++++++++++++++--
 src/equations/ideal_glm_mhd_1d.jl             |  8 -------
 src/equations/ideal_glm_mhd_2d.jl             |  9 --------
 src/equations/ideal_glm_mhd_3d.jl             | 10 ---------
 src/equations/polytropic_euler_2d.jl          |  9 --------
 src/equations/shallow_water_2d.jl             |  9 --------
 test/test_unit.jl                             |  2 --
 10 files changed, 20 insertions(+), 77 deletions(-)

diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
index 7884e65d6b..37a4f36210 100644
--- a/src/equations/compressible_euler_2d.jl
+++ b/src/equations/compressible_euler_2d.jl
@@ -1970,15 +1970,6 @@ end
     return v
 end
 
-@inline function velocity(u, normal_direction::AbstractVector,
-                          equations::CompressibleEulerEquations2D)
-    rho = u[1]
-    v1 = u[2] / rho
-    v2 = u[3] / rho
-    v = v1 * normal_direction[1] + v2 * normal_direction[2]
-    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 79ade4e7ae..56a24238fc 100644
--- a/src/equations/compressible_euler_3d.jl
+++ b/src/equations/compressible_euler_3d.jl
@@ -1706,16 +1706,6 @@ end
     return v
 end
 
-@inline function velocity(u, normal_direction::AbstractVector,
-                          equations::CompressibleEulerEquations3D)
-    rho = u[1]
-    v1 = u[2] / rho
-    v2 = u[3] / rho
-    v3 = u[4] / rho
-    v = v1 * normal_direction[1] + v2 * normal_direction[2] + v3 * normal_direction[3]
-    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_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl
index 65098f307b..2278fe58eb 100644
--- a/src/equations/compressible_euler_multicomponent_2d.jl
+++ b/src/equations/compressible_euler_multicomponent_2d.jl
@@ -839,13 +839,4 @@ end
     v = u[orientation] / rho
     return v
 end
-
-@inline function velocity(u, normal_direction::AbstractVector,
-                          equations::CompressibleEulerMulticomponentEquations2D)
-    rho = density(u, equations)
-    v1 = u[1] / rho
-    v2 = u[2] / rho
-    v = v1 * normal_direction[1] + v2 * normal_direction[2]
-    return v
-end
 end # @muladd
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index a95e619a3b..0304ea16cf 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -313,13 +313,31 @@ function prim2cons end
 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.
-
+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 24dd03a0d6..685ad6f8f3 100644
--- a/src/equations/ideal_glm_mhd_1d.jl
+++ b/src/equations/ideal_glm_mhd_1d.jl
@@ -569,14 +569,6 @@ end
     return v
 end
 
-@inline function velocity(u, normal_direction::AbstractVector,
-                          equations::IdealGlmMhdEquations1D)
-    rho = u[1]
-    v1 = u[2] / rho
-    v = v1 * normal_direction[1]
-    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 8d30920869..a2c748ae6f 100644
--- a/src/equations/ideal_glm_mhd_2d.jl
+++ b/src/equations/ideal_glm_mhd_2d.jl
@@ -1133,15 +1133,6 @@ end
     return v
 end
 
-@inline function velocity(u, normal_direction::AbstractVector,
-                          equations::IdealGlmMhdEquations2D)
-    rho = u[1]
-    v1 = u[2] / rho
-    v2 = u[3] / rho
-    v = v1 * normal_direction[1] + v2 * normal_direction[2]
-    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 5bffd23c84..1fbf67134c 100644
--- a/src/equations/ideal_glm_mhd_3d.jl
+++ b/src/equations/ideal_glm_mhd_3d.jl
@@ -1024,16 +1024,6 @@ end
     return v
 end
 
-@inline function velocity(u, normal_direction::AbstractVector,
-                          equations::IdealGlmMhdEquations3D)
-    rho = u[1]
-    v1 = u[2] / rho
-    v2 = u[3] / rho
-    v3 = u[4] / rho
-    v = v1 * normal_direction[1] + v2 * normal_direction[2] + v3 * normal_direction[3]
-    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 baa05de820..98f0f7f1dc 100644
--- a/src/equations/polytropic_euler_2d.jl
+++ b/src/equations/polytropic_euler_2d.jl
@@ -402,15 +402,6 @@ end
     return v
 end
 
-@inline function velocity(u, normal_direction::AbstractVector,
-                          equations::PolytropicEulerEquations2D)
-    rho = u[1]
-    v1 = u[2] / rho
-    v2 = u[3] / rho
-    v = v1 * normal_direction[1] + v2 * normal_direction[2]
-    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 3f5cc082dd..8b19d19a21 100644
--- a/src/equations/shallow_water_2d.jl
+++ b/src/equations/shallow_water_2d.jl
@@ -846,15 +846,6 @@ end
     return v
 end
 
-@inline function velocity(u, normal_direction::AbstractVector,
-                          equations::ShallowWaterEquations2D)
-    h = u[1]
-    v1 = u[2] / h
-    v2 = u[3] / h
-    v = v1 * normal_direction[1] + v2 * normal_direction[2]
-    return v
-end
-
 # Convert conservative variables to primitive
 @inline function cons2prim(u, equations::ShallowWaterEquations2D)
     h, _, _, b = u
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 7bf2448247..1cfe583795 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1811,8 +1811,6 @@ end
     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))
-    @test isapprox(velocity(u, SVector(normal_direction_2d[1]), equations_ideal_mhd_1d),
-                   v_normal_1d)
     for orientation in 1:3
         @test isapprox(velocity(u, orientation, equations_ideal_mhd_1d),
                        v_vector[orientation])