Implement max_abs_speed for real maximum speed

src/equations/shallow_water_multilayer_1d.jl

@@ -513,6 +513,27 @@ end
     return (max(abs(v_m_ll), abs(v_m_rr)) + max(c_ll, c_rr))
 end
 
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function Trixi.max_abs_speed(u_ll, u_rr,
+                                     orientation::Integer,
+                                     equations::ShallowWaterMultiLayerEquations1D)
+    # Unpack left and right state
+    h_ll = waterheight(u_ll, equations)
+    h_rr = waterheight(u_rr, equations)
+    hv_ll = momentum(u_ll, equations)
+    hv_rr = momentum(u_rr, equations)
+
+    # Get the averaged velocity
+    v_m_ll = sum(hv_ll) / sum(h_ll)
+    v_m_rr = sum(hv_rr) / sum(h_rr)
+
+    # Calculate the wave celerity on the left and right
+    c_ll = sqrt(equations.gravity * sum(h_ll))
+    c_rr = sqrt(equations.gravity * sum(h_rr))
+
+    return (max(abs(v_m_ll) + c_ll, abs(v_m_rr) + c_rr))
+end
+
 # Convert conservative variables to primitive
 @inline function Trixi.cons2prim(u, equations::ShallowWaterMultiLayerEquations1D)
     # Extract waterheight

src/equations/shallow_water_multilayer_2d.jl

@@ -796,6 +796,62 @@ end
             max(c_ll, c_rr) * norm(normal_direction))
 end
 
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function Trixi.max_abs_speed(u_ll, u_rr,
+                                     orientation::Integer,
+                                     equations::ShallowWaterMultiLayerEquations2D)
+    # Unpack left and right state
+    h_ll = waterheight(u_ll, equations)
+    h_rr = waterheight(u_rr, equations)
+
+    # Get the momentum quantities in the appropriate direction
+    if orientation == 1
+        h_v_ll, _ = momentum(u_ll, equations)
+        h_v_rr, _ = momentum(u_rr, equations)
+    else
+        _, h_v_ll = momentum(u_ll, equations)
+        _, h_v_rr = momentum(u_rr, equations)
+    end
+
+    # Get the averaged velocity
+    v_m_ll = sum(h_v_ll) / sum(h_ll)
+    v_m_rr = sum(h_v_rr) / sum(h_rr)
+
+    # Calculate the wave celerity on the left and right
+    c_ll = sqrt(equations.gravity * sum(h_ll))
+    c_rr = sqrt(equations.gravity * sum(h_rr))
+
+    return max(abs(v_m_ll) + c_ll, abs(v_m_rr) + c_rr)
+end
+
+@inline function Trixi.max_abs_speed(u_ll, u_rr,
+                                     normal_direction::AbstractVector,
+                                     equations::ShallowWaterMultiLayerEquations2D)
+    # Unpack left and right state
+    h_ll = waterheight(u_ll, equations)
+    h_rr = waterheight(u_rr, equations)
+    h_v1_ll, h_v2_ll = momentum(u_ll, equations)
+    h_v1_rr, h_v2_rr = momentum(u_rr, equations)
+
+    # Get the averaged velocity
+    v1_m_ll = sum(h_v1_ll) / sum(h_ll)
+    v2_m_ll = sum(h_v2_ll) / sum(h_ll)
+    v1_m_rr = sum(h_v1_rr) / sum(h_rr)
+    v2_m_rr = sum(h_v2_rr) / sum(h_rr)
+
+    # Compute velocity in the normal direction
+    v_dot_n_ll = v1_m_ll * normal_direction[1] + v2_m_ll * normal_direction[2]
+    v_dot_n_rr = v1_m_rr * normal_direction[1] + v2_m_rr * normal_direction[2]
+
+    # Calculate the wave celerity on the left and right
+    c_ll = sqrt(equations.gravity * sum(h_ll))
+    c_rr = sqrt(equations.gravity * sum(h_rr))
+
+    norm_ = norm(normal_direction)
+    # The normal velocities are already scaled by the norm
+    return (max(abs(v_dot_n_ll) + c_ll * norm_, abs(v_dot_n_rr)) + c_rr * norm_)
+end
+
 # Convert conservative variables to primitive
 @inline function Trixi.cons2prim(u, equations::ShallowWaterMultiLayerEquations2D)
     # Extract waterheight

src/equations/shallow_water_two_layer_1d.jl

@@ -393,6 +393,27 @@ end
     c_ll = sqrt(equations.gravity * (h_upper_ll + h_lower_ll))
     c_rr = sqrt(equations.gravity * (h_upper_rr + h_lower_rr))
 
+    return (max(abs(v_m_ll), abs(v_m_rr))) + max(c_ll, c_rr)
+end
+
+# Less "cautios", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function Trixi.max_abs_speed(u_ll, u_rr,
+                                     orientation::Integer,
+                                     equations::ShallowWaterTwoLayerEquations1D)
+    # Unpack left and right state
+    h_upper_ll, h_v_upper_ll, h_lower_ll, h_v_lower_ll, _ = u_ll
+    h_upper_rr, h_v_upper_rr, h_lower_rr, h_v_lower_rr, _ = u_rr
+
+    # Get the averaged velocity
+    v_m_ll = (h_v_upper_ll + h_v_lower_ll) / (h_upper_ll + h_lower_ll)
+    v_m_rr = (h_v_upper_rr + h_v_lower_rr) / (h_upper_rr + h_lower_rr)
+
+    # Calculate the wave celerity on the left and right
+    h_upper_ll, h_lower_ll = waterheight(u_ll, equations)
+    h_upper_rr, h_lower_rr = waterheight(u_rr, equations)
+    c_ll = sqrt(equations.gravity * (h_upper_ll + h_lower_ll))
+    c_rr = sqrt(equations.gravity * (h_upper_rr + h_lower_rr))
+
     return (max(abs(v_m_ll) + c_ll, abs(v_m_rr) + c_rr))
 end
 

src/equations/shallow_water_two_layer_2d.jl

@@ -669,6 +669,71 @@ end
     return max(abs(v_m_ll), abs(v_m_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 Trixi.max_abs_speed(u_ll, u_rr,
+                                     orientation::Integer,
+                                     equations::ShallowWaterTwoLayerEquations2D)
+    # Unpack left and right state
+    h_upper_ll, h_v1_upper_ll, h_v2_upper_ll, h_lower_ll, h_v1_lower_ll, h_v2_lower_ll, _ = u_ll
+    h_upper_rr, h_v1_upper_rr, h_v2_upper_rr, h_lower_rr, h_v1_lower_rr, h_v2_lower_rr, _ = u_rr
+
+    # Calculate averaged velocity of both layers
+    if orientation == 1
+        v_m_ll = (h_v1_upper_ll + h_v1_lower_ll) / (h_upper_ll + h_lower_ll)
+        v_m_rr = (h_v1_upper_rr + h_v1_lower_rr) / (h_upper_rr + h_lower_rr)
+    else
+        v_m_ll = (h_v2_upper_ll + h_v2_lower_ll) / (h_upper_ll + h_lower_ll)
+        v_m_rr = (h_v2_upper_rr + h_v2_lower_rr) / (h_upper_rr + h_lower_rr)
+    end
+
+    # Calculate the wave celerity on the left and right
+    h_upper_ll, h_lower_ll = waterheight(u_ll, equations)
+    h_upper_rr, h_lower_rr = waterheight(u_rr, equations)
+
+    c_ll = sqrt(equations.gravity * (h_upper_ll + h_lower_ll))
+    c_rr = sqrt(equations.gravity * (h_upper_rr + h_lower_rr))
+
+    return max(abs(v_m_ll) + c_ll, abs(v_m_rr) + c_rr)
+end
+
+@inline function Trixi.max_abs_speed(u_ll, u_rr,
+                                     normal_direction::AbstractVector,
+                                     equations::ShallowWaterTwoLayerEquations2D)
+    # Unpack left and right state
+    h_upper_ll, _, _, h_lower_ll, _, _, _ = u_ll
+    h_upper_rr, _, _, h_lower_rr, _, _, _ = u_rr
+
+    # Extract and compute the velocities in the normal direction
+    v1_upper_ll, v2_upper_ll, v1_lower_ll, v2_lower_ll = velocity(u_ll, equations)
+    v1_upper_rr, v2_upper_rr, v1_lower_rr, v2_lower_rr = velocity(u_rr, equations)
+
+    v_upper_dot_n_ll = v1_upper_ll * normal_direction[1] +
+                       v2_upper_ll * normal_direction[2]
+    v_upper_dot_n_rr = v1_upper_rr * normal_direction[1] +
+                       v2_upper_rr * normal_direction[2]
+    v_lower_dot_n_ll = v1_lower_ll * normal_direction[1] +
+                       v2_lower_ll * normal_direction[2]
+    v_lower_dot_n_rr = v1_lower_rr * normal_direction[1] +
+                       v2_lower_rr * normal_direction[2]
+
+    # Calculate averaged velocity of both layers
+    v_m_ll = (v_upper_dot_n_ll * h_upper_ll + v_lower_dot_n_ll * h_lower_ll) /
+             (h_upper_ll + h_lower_ll)
+    v_m_rr = (v_upper_dot_n_rr * h_upper_rr + v_lower_dot_n_rr * h_lower_rr) /
+             (h_upper_rr + h_lower_rr)
+
+    # Compute the wave celerity on the left and right
+    h_upper_ll, h_lower_ll = waterheight(u_ll, equations)
+    h_upper_rr, h_lower_rr = waterheight(u_rr, equations)
+
+    c_ll = sqrt(equations.gravity * (h_upper_ll + h_lower_ll))
+    c_rr = sqrt(equations.gravity * (h_upper_rr + h_lower_rr))
+
+    norm_ = norm(normal_direction)
+    # The normal velocities are already scaled by the norm
+    return (max(abs(v_m_ll) + c_ll * norm_, abs(v_m_rr) + c_rr * norm_))
+end
+
 # Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography
 @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr,
                                                               orientation_or_normal_direction,