From ba296e680eed5ac6de4ec9a8151505d9e46094c7 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 20 Jan 2025 13:32:21 +0100 Subject: [PATCH] Implement max_abs_speed for real maximum speed --- src/equations/shallow_water_multilayer_1d.jl | 21 +++++++ src/equations/shallow_water_multilayer_2d.jl | 56 +++++++++++++++++ src/equations/shallow_water_two_layer_1d.jl | 21 +++++++ src/equations/shallow_water_two_layer_2d.jl | 65 ++++++++++++++++++++ 4 files changed, 163 insertions(+) diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl index d9e18b2..753fde2 100644 --- a/src/equations/shallow_water_multilayer_1d.jl +++ b/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 diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl index ad59766..a226484 100644 --- a/src/equations/shallow_water_multilayer_2d.jl +++ b/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 diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index a3a0a83..89359fa 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/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 diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index 12abe85..db19ce6 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/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,