Implement max_abs_speed for real maximum speed

src/Trixi.jl

@@ -190,7 +190,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
        hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal,
        FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs,
        DissipationLaxFriedrichsEntropyVariables,
-       FluxLaxFriedrichs, max_abs_speed_naive,
+       FluxLaxFriedrichs, max_abs_speed_naive, max_abs_speed,
        FluxHLL, min_max_speed_naive, min_max_speed_davis, min_max_speed_einfeldt,
        FluxLMARS,
        FluxRotated,

src/equations/acoustic_perturbation_2d.jl

@@ -298,6 +298,24 @@ end
     λ_max = max(abs(v_ll), abs(v_rr)) + max(c_mean_ll, c_mean_rr)
 end
 
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
+                               equations::AcousticPerturbationEquations2D)
+    # Calculate v = v_prime + v_mean
+    v_prime_ll = u_ll[orientation]
+    v_prime_rr = u_rr[orientation]
+    v_mean_ll = u_ll[orientation + 3]
+    v_mean_rr = u_rr[orientation + 3]
+
+    v_ll = v_prime_ll + v_mean_ll
+    v_rr = v_prime_rr + v_mean_rr
+
+    c_mean_ll = u_ll[6]
+    c_mean_rr = u_rr[6]
+
+    λ_max = max(abs(v_ll) + c_mean_ll, abs(v_rr) + c_mean_rr)
+end
+
 # Calculate 1D flux for a single point in the normal direction
 # Note, this directional vector is not normalized
 @inline function flux(u, normal_direction::AbstractVector,
@@ -341,6 +359,26 @@ end
             max(c_mean_ll, c_mean_rr) * norm(normal_direction)
 end
 
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, normal_direction::AbstractVector,
+                               equations::AcousticPerturbationEquations2D)
+    # Calculate v = v_prime + v_mean
+    v_prime_ll = normal_direction[1] * u_ll[1] + normal_direction[2] * u_ll[2]
+    v_prime_rr = normal_direction[1] * u_rr[1] + normal_direction[2] * u_rr[2]
+    v_mean_ll = normal_direction[1] * u_ll[4] + normal_direction[2] * u_ll[5]
+    v_mean_rr = normal_direction[1] * u_rr[4] + normal_direction[2] * u_rr[5]
+
+    v_ll = v_prime_ll + v_mean_ll
+    v_rr = v_prime_rr + v_mean_rr
+
+    c_mean_ll = u_ll[6]
+    c_mean_rr = u_rr[6]
+
+    norm_ = norm(normal_direction)
+    # The v_normals are already scaled by the norm
+    λ_max = max(abs(v_ll) + c_mean_ll * norm_, abs(v_rr) + c_mean_rr * norm_)
+end
+
 # Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the mean values
 @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr,
                                                               orientation_or_normal_direction,

src/equations/compressible_euler_1d.jl

@@ -720,6 +720,25 @@ end
     return λ_min, λ_max
 end
 
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
+                               equations::CompressibleEulerEquations1D)
+    rho_ll, rho_v1_ll, rho_e_ll = u_ll
+    rho_rr, rho_v1_rr, rho_e_rr = u_rr
+
+    # Calculate primitive variables and speed of sound
+    v1_ll = rho_v1_ll / rho_ll
+    v_mag_ll = abs(v1_ll)
+    p_ll = (equations.gamma - 1) * (rho_e_ll - 0.5f0 * rho_ll * v_mag_ll^2)
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    v1_rr = rho_v1_rr / rho_rr
+    v_mag_rr = abs(v1_rr)
+    p_rr = (equations.gamma - 1) * (rho_e_rr - 0.5f0 * rho_rr * v_mag_rr^2)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    λ_max = max(v_mag_ll + c_ll, v_mag_rr + c_rr)
+end
+
 # More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
 @inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
                                      equations::CompressibleEulerEquations1D)

src/equations/compressible_euler_2d.jl

@@ -1433,6 +1433,50 @@ end
     return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction)
 end
 
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
+                               equations::CompressibleEulerEquations2D)
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+
+    # Get the velocity value in the appropriate direction
+    if orientation == 1
+        v_ll = v1_ll
+        v_rr = v1_rr
+    else # orientation == 2
+        v_ll = v2_ll
+        v_rr = v2_rr
+    end
+    # Calculate sound speeds
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    λ_max = max(abs(v_ll) + c_ll, abs(v_rr) + c_rr)
+end
+
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, normal_direction::AbstractVector,
+                               equations::CompressibleEulerEquations2D)
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+
+    # Calculate normal velocities and sound speed
+    # left
+    v_ll = (v1_ll * normal_direction[1]
+            +
+            v2_ll * normal_direction[2])
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    # right
+    v_rr = (v1_rr * normal_direction[1]
+            +
+            v2_rr * normal_direction[2])
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    norm_ = norm(normal_direction)
+    return max(abs(v_ll) + c_ll * norm_,
+               abs(v_rr) + c_rr * norm_)
+end
+
 # Calculate estimate for minimum and maximum wave speeds for HLL-type fluxes
 @inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::CompressibleEulerEquations2D)

src/equations/compressible_euler_3d.jl

@@ -1133,6 +1133,52 @@ end
     return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction)
 end
 
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
+                               equations::CompressibleEulerEquations3D)
+    rho_ll, v1_ll, v2_ll, v3_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, v3_rr, p_rr = cons2prim(u_rr, equations)
+
+    # Get the velocity value in the appropriate direction
+    if orientation == 1
+        v_ll = v1_ll
+        v_rr = v1_rr
+    elseif orientation == 2
+        v_ll = v2_ll
+        v_rr = v2_rr
+    else # orientation == 3
+        v_ll = v3_ll
+        v_rr = v3_rr
+    end
+    # Calculate sound speeds
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    λ_max = max(abs(v_ll) + c_ll, abs(v_rr) + c_rr)
+end
+
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, normal_direction::AbstractVector,
+                               equations::CompressibleEulerEquations3D)
+    rho_ll, v1_ll, v2_ll, v3_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, v3_rr, p_rr = cons2prim(u_rr, equations)
+
+    # Calculate normal velocities and sound speed
+    # left
+    v_ll = (v1_ll * normal_direction[1]
+            + v2_ll * normal_direction[2]
+            + v3_ll * normal_direction[3])
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    # right
+    v_rr = (v1_rr * normal_direction[1]
+            + v2_rr * normal_direction[2]
+            + v3_rr * normal_direction[3])
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    norm_ = norm(normal_direction)
+    return max(abs(v_ll) + c_ll * norm_, abs(v_rr) + c_rr * norm_)
+end
+
 # Calculate estimates for minimum and maximum wave speeds for HLL-type fluxes
 @inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::CompressibleEulerEquations3D)

src/equations/compressible_euler_multicomponent_1d.jl

@@ -408,6 +408,29 @@ end
     λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr)
 end
 
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
+                               equations::CompressibleEulerMulticomponentEquations1D)
+    rho_v1_ll, rho_e_ll = u_ll
+    rho_v1_rr, rho_e_rr = u_rr
+
+    # Calculate primitive variables and speed of sound
+    rho_ll = density(u_ll, equations)
+    rho_rr = density(u_rr, equations)
+    gamma_ll = totalgamma(u_ll, equations)
+    gamma_rr = totalgamma(u_rr, equations)
+
+    v_ll = rho_v1_ll / rho_ll
+    v_rr = rho_v1_rr / rho_rr
+
+    p_ll = (gamma_ll - 1) * (rho_e_ll - 0.5f0 * rho_ll * v_ll^2)
+    p_rr = (gamma_rr - 1) * (rho_e_rr - 0.5f0 * rho_rr * v_rr^2)
+    c_ll = sqrt(gamma_ll * p_ll / rho_ll)
+    c_rr = sqrt(gamma_rr * p_rr / rho_rr)
+
+    λ_max = max(abs(v_ll) + c_ll, abs(v_rr) + c_rr)
+end
+
 @inline function max_abs_speeds(u,
                                 equations::CompressibleEulerMulticomponentEquations1D)
     rho_v1, rho_e = u

src/equations/compressible_euler_multicomponent_2d.jl

@@ -567,6 +567,36 @@ end
     λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr)
 end
 
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
+                               equations::CompressibleEulerMulticomponentEquations2D)
+    rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
+    rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr
+
+    # Get the density and gas gamma
+    rho_ll = density(u_ll, equations)
+    rho_rr = density(u_rr, equations)
+    gamma_ll = totalgamma(u_ll, equations)
+    gamma_rr = totalgamma(u_rr, equations)
+
+    # Get the velocities based on direction
+    if orientation == 1
+        v_ll = rho_v1_ll / rho_ll
+        v_rr = rho_v1_rr / rho_rr
+    else # orientation == 2
+        v_ll = rho_v2_ll / rho_ll
+        v_rr = rho_v2_rr / rho_rr
+    end
+
+    # Compute the sound speeds on the left and right
+    p_ll = (gamma_ll - 1) * (rho_e_ll - 0.5f0 * (rho_v1_ll^2 + rho_v2_ll^2) / rho_ll)
+    c_ll = sqrt(gamma_ll * p_ll / rho_ll)
+    p_rr = (gamma_rr - 1) * (rho_e_rr - 0.5f0 * (rho_v1_rr^2 + rho_v2_rr^2) / rho_rr)
+    c_rr = sqrt(gamma_rr * p_rr / rho_rr)
+
+    λ_max = max(abs(v_ll) + c_ll, abs(v_rr) + c_rr)
+end
+
 @inline function max_abs_speeds(u,
                                 equations::CompressibleEulerMulticomponentEquations2D)
     rho_v1, rho_v2, rho_e = u

src/equations/compressible_euler_quasi_1d.jl

@@ -288,6 +288,29 @@ end
     λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
 end
 
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
+                               equations::CompressibleEulerEquationsQuasi1D)
+    a_rho_ll, a_rho_v1_ll, a_e_ll, a_ll = u_ll
+    a_rho_rr, a_rho_v1_rr, a_e_rr, a_rr = u_rr
+
+    # Calculate primitive variables and speed of sound
+    rho_ll = a_rho_ll / a_ll
+    e_ll = a_e_ll / a_ll
+    v1_ll = a_rho_v1_ll / a_rho_ll
+    v_mag_ll = abs(v1_ll)
+    p_ll = (equations.gamma - 1) * (e_ll - 0.5f0 * rho_ll * v_mag_ll^2)
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    rho_rr = a_rho_rr / a_rr
+    e_rr = a_e_rr / a_rr
+    v1_rr = a_rho_v1_rr / a_rho_rr
+    v_mag_rr = abs(v1_rr)
+    p_rr = (equations.gamma - 1) * (e_rr - 0.5f0 * rho_rr * v_mag_rr^2)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    λ_max = max(v_mag_ll + c_ll, v_mag_rr + c_rr)
+end
+
 @inline function max_abs_speeds(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
     rho = a_rho / a

src/equations/ideal_glm_mhd_1d.jl

@@ -426,6 +426,24 @@ end
     λ_max = max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr)
 end
 
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
+                               equations::IdealGlmMhdEquations1D)
+    rho_ll, rho_v1_ll, _ = u_ll
+    rho_rr, rho_v1_rr, _ = u_rr
+
+    # Calculate velocities (ignore orientation since it is always "1" in 1D)
+    # and fast magnetoacoustic wave speeds
+    # left
+    v_ll = rho_v1_ll / rho_ll
+    cf_ll = calc_fast_wavespeed(u_ll, orientation, equations)
+    # right
+    v_rr = rho_v1_rr / rho_rr
+    cf_rr = calc_fast_wavespeed(u_rr, orientation, equations)
+
+    λ_max = max(abs(v_ll) + cf_ll, abs(v_rr) + cf_rr)
+end
+
 # Calculate estimates for minimum and maximum wave speeds for HLL-type fluxes
 @inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::IdealGlmMhdEquations1D)

src/equations/ideal_glm_mhd_2d.jl

@@ -780,6 +780,53 @@ end
     return max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr)
 end
 
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
+                               equations::IdealGlmMhdEquations2D)
+    rho_ll, rho_v1_ll, rho_v2_ll, _ = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, _ = u_rr
+
+    # Calculate the left/right velocities and fast magnetoacoustic wave speeds
+    if orientation == 1
+        v_ll = rho_v1_ll / rho_ll
+        v_rr = rho_v1_rr / rho_rr
+    else # orientation == 2
+        v_ll = rho_v2_ll / rho_ll
+        v_rr = rho_v2_rr / rho_rr
+    end
+    cf_ll = calc_fast_wavespeed(u_ll, orientation, equations)
+    cf_rr = calc_fast_wavespeed(u_rr, orientation, equations)
+
+    return max(abs(v_ll) + cf_ll, abs(v_rr) + cf_rr)
+end
+
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, normal_direction::AbstractVector,
+                               equations::IdealGlmMhdEquations2D)
+    # return max(v_mag_ll, v_mag_rr) + max(cf_ll, cf_rr)
+    rho_ll, rho_v1_ll, rho_v2_ll, _ = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, _ = u_rr
+
+    # Calculate normal velocities and fast magnetoacoustic wave speeds
+    # left
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+    v_ll = (v1_ll * normal_direction[1]
+            +
+            v2_ll * normal_direction[2])
+    cf_ll = calc_fast_wavespeed(u_ll, normal_direction, equations)
+    # right
+    v1_rr = rho_v1_rr / rho_rr
+    v2_rr = rho_v2_rr / rho_rr
+    v_rr = (v1_rr * normal_direction[1]
+            +
+            v2_rr * normal_direction[2])
+    cf_rr = calc_fast_wavespeed(u_rr, normal_direction, equations)
+
+    # wave speeds already scaled by norm(normal_direction) in [`calc_fast_wavespeed`](@ref)
+    return max(abs(v_ll) + cf_ll, abs(v_rr) + cf_rr)
+end
+
 # Calculate estimate for minimum and maximum wave speeds for HLL-type fluxes
 @inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::IdealGlmMhdEquations2D)