Add velocity functions for different equations

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,

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)

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)

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)

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

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

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)
 

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,30 @@ 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 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)
     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

src/equations/ideal_glm_mhd_2d.jl

@@ -1119,6 +1119,30 @@ 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 velocity(u, normal_direction::AbstractVector,
+                          equations::IdealGlmMhdEquations2D)
+    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::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

src/equations/ideal_glm_mhd_3d.jl

@@ -1010,6 +1010,30 @@ 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

src/equations/polytropic_euler_2d.jl

@@ -389,6 +389,28 @@ 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 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

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,