From 03eb72530d86680bd9cf43cdba0b0a44176a75db Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 21 Dec 2023 15:37:38 +0100 Subject: [PATCH] add HLLC flux for non-cartesian meshes --- src/equations/compressible_euler_2d.jl | 101 ++++++++++++++++++++- src/equations/compressible_euler_3d.jl | 117 +++++++++++++++++++++++-- 2 files changed, 206 insertions(+), 12 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index a992f99eaf4..8f369f82a5e 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1185,18 +1185,18 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, if orientation == 1 # x-direction vel_L = v1_ll vel_R = v1_rr - ekin_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2 elseif orientation == 2 # y-direction vel_L = v2_ll vel_R = v2_rr - ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2 end vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho - ekin_roe = 0.5 * (vel_roe^2 + ekin_roe / sum_sqrt_rho^2) + v1_roe = sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr + v2_roe = sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr + vel_roe_mag = (v1_roe^2 + v2_roe^2) / sum_sqrt_rho^2 H_ll = (rho_e_ll + p_ll) / rho_ll H_rr = (rho_e_rr + p_rr) / rho_rr H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - ekin_roe)) + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) Ssl = min(vel_L - c_ll, vel_roe - c_roe) Ssr = max(vel_R + c_rr, vel_roe + c_roe) sMu_L = Ssl - vel_L @@ -1252,6 +1252,99 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, return SVector(f1, f2, f3, f4) end +function flux_hllc(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + # Calculate primitive variables and speed of sound + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + norm_ = norm(normal_direction) + norm_sq = norm_ * norm_ + inv_norm_sq = 1.0 / norm_sq + + c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_ + c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_ + + # Obtain left and right fluxes + f_ll = flux(u_ll, normal_direction, equations) + f_rr = flux(u_rr, normal_direction, equations) + + # Compute Roe averages + sqrt_rho_ll = sqrt(rho_ll) + sqrt_rho_rr = sqrt(rho_rr) + sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr + + v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) / sum_sqrt_rho + v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) / sum_sqrt_rho + vel_roe = v1_roe * normal_direction[1] + v2_roe * normal_direction[2] + vel_roe_mag = v1_roe^2 + v2_roe^2 + + e_ll = u_ll[4] / rho_ll + e_rr = u_rr[4] / rho_rr + + H_ll = (u_ll[4] + p_ll) / rho_ll + H_rr = (u_rr[4] + p_rr) / rho_rr + + H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) * norm_ + + Ssl = min(v_dot_n_ll - c_ll, vel_roe - c_roe) + Ssr = max(v_dot_n_rr + c_rr, vel_roe + c_roe) + sMu_L = Ssl - v_dot_n_ll + sMu_R = Ssr - v_dot_n_rr + + if Ssl >= 0.0 + f1 = f_ll[1] + f2 = f_ll[2] + f3 = f_ll[3] + f4 = f_ll[4] + elseif Ssr <= 0.0 + f1 = f_rr[1] + f2 = f_rr[2] + f3 = f_rr[3] + f4 = f_rr[4] + else + SStar = (rho_ll * v_dot_n_ll * sMu_L - rho_rr * v_dot_n_rr * sMu_R + + + (p_rr - p_ll) * norm_sq) / (rho_ll * sMu_L - rho_rr * sMu_R) + if Ssl <= 0.0 <= SStar + densStar = rho_ll * sMu_L / (Ssl - SStar) + enerStar = e_ll + + (SStar - v_dot_n_ll) * + (SStar * inv_norm_sq + p_ll / (rho_ll * sMu_L)) + UStar1 = densStar + UStar2 = densStar * + (v1_ll + (SStar - v_dot_n_ll) * normal_direction[1] * inv_norm_sq) + UStar3 = densStar * + (v2_ll + (SStar - v_dot_n_ll) * normal_direction[2] * inv_norm_sq) + UStar4 = densStar * enerStar + f1 = f_ll[1] + Ssl * (UStar1 - u_ll[1]) + f2 = f_ll[2] + Ssl * (UStar2 - u_ll[2]) + f3 = f_ll[3] + Ssl * (UStar3 - u_ll[3]) + f4 = f_ll[4] + Ssl * (UStar4 - u_ll[4]) + else + densStar = rho_rr * sMu_R / (Ssr - SStar) + enerStar = e_rr + + (SStar - v_dot_n_rr) * + (SStar * inv_norm_sq + p_rr / (rho_rr * sMu_R)) + UStar1 = densStar + UStar2 = densStar * + (v1_rr + (SStar - v_dot_n_rr) * normal_direction[1] * inv_norm_sq) + UStar3 = densStar * + (v2_rr + (SStar - v_dot_n_rr) * normal_direction[2] * inv_norm_sq) + UStar4 = densStar * enerStar + f1 = f_rr[1] + Ssr * (UStar1 - u_rr[1]) + f2 = f_rr[2] + Ssr * (UStar2 - u_rr[2]) + f3 = f_rr[3] + Ssr * (UStar3 - u_rr[3]) + f4 = f_rr[4] + Ssr * (UStar4 - u_rr[4]) + end + end + return SVector(f1, f2, f3, f4) +end + """ min_max_speed_einfeldt(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index fc56f58025b..d6de349fd91 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -1231,25 +1231,22 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, if orientation == 1 # x-direction vel_L = v1_ll vel_R = v1_rr - ekin_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2 + - (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr)^2 elseif orientation == 2 # y-direction vel_L = v2_ll vel_R = v2_rr - ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2 + - (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr)^2 else # z-direction vel_L = v3_ll vel_R = v3_rr - ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2 + - (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2 end vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho - ekin_roe = 0.5 * (vel_roe^2 + ekin_roe / sum_sqrt_rho^2) + v1_roe = sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr + v2_roe = sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr + v3_roe = sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr + vel_roe_mag = (v1_roe^2 + v2_roe^2 + v3_roe^2) / sum_sqrt_rho^2 H_ll = (rho_e_ll + p_ll) / rho_ll H_rr = (rho_e_rr + p_rr) / rho_rr H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - ekin_roe)) + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) Ssl = min(vel_L - c_ll, vel_roe - c_roe) Ssr = max(vel_R + c_rr, vel_roe + c_roe) sMu_L = Ssl - vel_L @@ -1321,6 +1318,110 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, return SVector(f1, f2, f3, f4, f5) end +function flux_hllc(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations3D) + # Calculate primitive variables and speed of sound + 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) + + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3] + + norm_ = norm(normal_direction) + norm_sq = norm_ * norm_ + inv_norm_sq = 1.0 / norm_sq + + c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_ + c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_ + + # Obtain left and right fluxes + f_ll = flux(u_ll, normal_direction, equations) + f_rr = flux(u_rr, normal_direction, equations) + + # Compute Roe averages + sqrt_rho_ll = sqrt(rho_ll) + sqrt_rho_rr = sqrt(rho_rr) + sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr + + v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) / sum_sqrt_rho + v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) / sum_sqrt_rho + v3_roe = (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr) / sum_sqrt_rho + vel_roe = v1_roe * normal_direction[1] + v2_roe * normal_direction[2] + + v3_roe * normal_direction[3] + vel_roe_mag = v1_roe^2 + v2_roe^2 + v3_roe^2 + + e_ll = u_ll[5] / rho_ll + e_rr = u_rr[5] / rho_rr + + H_ll = (u_ll[5] + p_ll) / rho_ll + H_rr = (u_rr[5] + p_rr) / rho_rr + + H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) * norm_ + + Ssl = min(v_dot_n_ll - c_ll, vel_roe - c_roe) + Ssr = max(v_dot_n_rr + c_rr, vel_roe + c_roe) + sMu_L = Ssl - v_dot_n_ll + sMu_R = Ssr - v_dot_n_rr + + if Ssl >= 0.0 + f1 = f_ll[1] + f2 = f_ll[2] + f3 = f_ll[3] + f4 = f_ll[4] + f5 = f_ll[5] + elseif Ssr <= 0.0 + f1 = f_rr[1] + f2 = f_rr[2] + f3 = f_rr[3] + f4 = f_rr[4] + f5 = f_rr[5] + else + SStar = (rho_ll * v_dot_n_ll * sMu_L - rho_rr * v_dot_n_rr * sMu_R + + (p_rr - p_ll) * norm_sq) / (rho_ll * sMu_L - rho_rr * sMu_R) + if Ssl <= 0.0 <= SStar + densStar = rho_ll * sMu_L / (Ssl - SStar) + enerStar = e_ll + + (SStar - v_dot_n_ll) * + (SStar * inv_norm_sq + p_ll / (rho_ll * sMu_L)) + UStar1 = densStar + UStar2 = densStar * + (v1_ll + (SStar - v_dot_n_ll) * normal_direction[1] * inv_norm_sq) + UStar3 = densStar * + (v2_ll + (SStar - v_dot_n_ll) * normal_direction[2] * inv_norm_sq) + UStar4 = densStar * + (v3_ll + (SStar - v_dot_n_ll) * normal_direction[3] * inv_norm_sq) + UStar5 = densStar * enerStar + f1 = f_ll[1] + Ssl * (UStar1 - u_ll[1]) + f2 = f_ll[2] + Ssl * (UStar2 - u_ll[2]) + f3 = f_ll[3] + Ssl * (UStar3 - u_ll[3]) + f4 = f_ll[4] + Ssl * (UStar4 - u_ll[4]) + f5 = f_ll[5] + Ssl * (UStar5 - u_ll[5]) + else + densStar = rho_rr * sMu_R / (Ssr - SStar) + enerStar = e_rr + + (SStar - v_dot_n_rr) * + (SStar * inv_norm_sq + p_rr / (rho_rr * sMu_R)) + UStar1 = densStar + UStar2 = densStar * + (v1_rr + (SStar - v_dot_n_rr) * normal_direction[1] * inv_norm_sq) + UStar3 = densStar * + (v2_rr + (SStar - v_dot_n_rr) * normal_direction[2] * inv_norm_sq) + UStar4 = densStar * + (v3_rr + (SStar - v_dot_n_rr) * normal_direction[3] * inv_norm_sq) + UStar5 = densStar * enerStar + f1 = f_rr[1] + Ssr * (UStar1 - u_rr[1]) + f2 = f_rr[2] + Ssr * (UStar2 - u_rr[2]) + f3 = f_rr[3] + Ssr * (UStar3 - u_rr[3]) + f4 = f_rr[4] + Ssr * (UStar4 - u_rr[4]) + f5 = f_rr[5] + Ssr * (UStar5 - u_rr[5]) + end + end + return SVector(f1, f2, f3, f4, f5) +end + """ min_max_speed_einfeldt(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D)