From fdbd01261e1d9ef1779eed6b752784585f53e9a1 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 20 Jun 2024 22:02:50 -1000 Subject: [PATCH 1/6] start --- .../ideal_glm_mhd_multicomponent_1d.jl | 93 +++++------ .../ideal_glm_mhd_multicomponent_2d.jl | 154 +++++++++--------- src/equations/inviscid_burgers_1d.jl | 18 +- 3 files changed, 134 insertions(+), 131 deletions(-) diff --git a/src/equations/ideal_glm_mhd_multicomponent_1d.jl b/src/equations/ideal_glm_mhd_multicomponent_1d.jl index dad7c27e86c..7ed3fef9d82 100644 --- a/src/equations/ideal_glm_mhd_multicomponent_1d.jl +++ b/src/equations/ideal_glm_mhd_multicomponent_1d.jl @@ -120,7 +120,7 @@ function initial_condition_weak_blast_wave(x, t, phi = atan(x_norm) # Calculate primitive variables - if r > 0.5 + if r > 0.5f0 rho = 1.0 prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / (1 - @@ -135,8 +135,8 @@ function initial_condition_weak_blast_wave(x, t, rho for i in eachcomponent(equations)) end - v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi) - p = r > 0.5 ? 1.0 : 1.245 + v1 = r > 0.5f0 ? 0.0 : 0.1882 * cos(phi) + p = r > 0.5f0 ? 1.0 : 1.245 prim_other = SVector{7, real(equations)}(v1, 0.0, 0.0, p, 1.0, 1.0, 1.0) @@ -153,8 +153,8 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * rho * (v1^2 + v2^2 + v3^2) - mag_en = 0.5 * (B1^2 + B2^2 + B3^2) + kin_en = 0.5f0 * rho * (v1^2 + v2^2 + v3^2) + mag_en = 0.5f0 * (B1^2 + B2^2 + B3^2) gamma = totalgamma(u, equations) p = (gamma - 1) * (rho_e - kin_en - mag_en) @@ -198,7 +198,7 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 7], u_rr[i + 7]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 7] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 7] + u_rr[i + 7]) for i in eachcomponent(equations)) @@ -217,16 +217,16 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, vel_dot_mag_rr = v1_rr * B1_rr + v2_rr * B2_rr + v3_rr * B3_rr # Compute the necessary mean values needed for either direction - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) v_sum = v1_avg + v2_avg + v3_avg - B1_avg = 0.5 * (B1_ll + B1_rr) - B2_avg = 0.5 * (B2_ll + B2_rr) - B3_avg = 0.5 * (B3_ll + B3_rr) - vel_norm_avg = 0.5 * (vel_norm_ll + vel_norm_rr) - mag_norm_avg = 0.5 * (mag_norm_ll + mag_norm_rr) - vel_dot_mag_avg = 0.5 * (vel_dot_mag_ll + vel_dot_mag_rr) + B1_avg = 0.5f0 * (B1_ll + B1_rr) + B2_avg = 0.5f0 * (B2_ll + B2_rr) + B3_avg = 0.5f0 * (B3_ll + B3_rr) + vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr) + mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr) + vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr) enth = zero(v_sum) help1_ll = zero(v1_ll) @@ -238,9 +238,9 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, help1_rr += u_rr[i + 7] * cv[i] end - T_ll = (rho_e_ll - 0.5 * rho_ll * (vel_norm_ll) - 0.5 * mag_norm_ll) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (vel_norm_rr) - 0.5 * mag_norm_rr) / help1_rr - T = 0.5 * (1.0 / T_ll + 1.0 / T_rr) + T_ll = (rho_e_ll - 0.5f0 * rho_ll * (vel_norm_ll) - 0.5f0 * mag_norm_ll) / help1_ll + T_rr = (rho_e_rr - 0.5f0 * rho_rr * (vel_norm_rr) - 0.5f0 * mag_norm_rr) / help1_rr + T = 0.5f0 * (1.0 / T_ll + 1.0 / T_rr) T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) # Calculate fluxes depending on orientation with specific direction averages @@ -253,7 +253,7 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, help1 += f_rho[i] * cv[i] help2 += f_rho[i] end - f1 = help2 * v1_avg + enth / T + 0.5 * mag_norm_avg - B1_avg * B1_avg + f1 = help2 * v1_avg + enth / T + 0.5f0 * mag_norm_avg - B1_avg * B1_avg f2 = help2 * v2_avg - B1_avg * B2_avg f3 = help2 * v3_avg - B1_avg * B3_avg f5 = 0.0 @@ -261,11 +261,12 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, f7 = v1_avg * B3_avg - v3_avg * B1_avg # total energy flux is complicated and involves the previous eight components - v1_mag_avg = 0.5 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) + v1_mag_avg = 0.5f0 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) - f4 = (help1 / T_log) - 0.5 * (vel_norm_avg) * (help2) + f1 * v1_avg + f2 * v2_avg + + f4 = (help1 / T_log) - 0.5f0 * (vel_norm_avg) * (help2) + f1 * v1_avg + + f2 * v2_avg + f3 * v3_avg + - f5 * B1_avg + f6 * B2_avg + f7 * B3_avg - 0.5 * v1_mag_avg + + f5 * B1_avg + f6 * B2_avg + f7 * B3_avg - 0.5f0 * v1_mag_avg + B1_avg * vel_dot_mag_avg f_other = SVector{7, real(equations)}(f1, f2, f3, f4, f5, f6, f7) @@ -309,19 +310,19 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations. # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_avg = 0.5 * (p_ll + p_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) - magnetic_square_avg = 0.5 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) - inv_gamma_minus_one = 1 / (totalgamma(0.5 * (u_ll + u_rr), equations) - 1) + inv_gamma_minus_one = 1 / (totalgamma(0.5f0 * (u_ll + u_rr), equations) - 1) rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 7], u_rr[i + 7]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 7] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 7] + u_rr[i + 7]) for i in eachcomponent(equations)) @@ -334,17 +335,17 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations. # Calculate fluxes depending on orientation with specific direction averages f2 = f1 * v1_avg + p_avg + magnetic_square_avg - - 0.5 * (B1_ll * B1_rr + B1_rr * B1_ll) - f3 = f1 * v2_avg - 0.5 * (B1_ll * B2_rr + B1_rr * B2_ll) - f4 = f1 * v3_avg - 0.5 * (B1_ll * B3_rr + B1_rr * B3_ll) + 0.5f0 * (B1_ll * B1_rr + B1_rr * B1_ll) + f3 = f1 * v2_avg - 0.5f0 * (B1_ll * B2_rr + B1_rr * B2_ll) + f4 = f1 * v3_avg - 0.5f0 * (B1_ll * B3_rr + B1_rr * B3_ll) #f5 below f6 = 0.0 - f7 = 0.5 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) - f8 = 0.5 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) + f7 = 0.5f0 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) + f8 = 0.5f0 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (+p_ll * v1_rr + p_rr * v1_ll + 0.5f0 * (+p_ll * v1_rr + p_rr * v1_ll + (v1_ll * B2_ll * B2_rr + v1_rr * B2_rr * B2_ll) + (v1_ll * B3_ll * B3_rr + v1_rr * B3_rr * B3_ll) - @@ -405,7 +406,7 @@ function cons2prim(u, equations::IdealGlmMhdMulticomponentEquations1D) gamma = totalgamma(u, equations) p = (gamma - 1) * - (rho_e - 0.5 * rho * (v1^2 + v2^2 + v3^2) - 0.5 * (B1^2 + B2^2 + B3^2)) + (rho_e - 0.5f0 * rho * (v1^2 + v2^2 + v3^2) - 0.5f0 * (B1^2 + B2^2 + B3^2)) prim_other = SVector{7, real(equations)}(v1, v2, v3, p, B1, B2, B3) return vcat(prim_other, prim_rho) end @@ -422,7 +423,7 @@ end v3 = rho_v3 / rho v_square = v1^2 + v2^2 + v3^2 gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * v_square - 0.5 * (B1^2 + B2^2 + B3^2)) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2)) s = log(p) - gamma * log(rho) rho_p = rho / p @@ -433,7 +434,7 @@ end help1 += u[i + 7] * cv[i] end - T = (rho_e - 0.5 * rho * v_square - 0.5 * (B1^2 + B2^2 + B3^2)) / (help1) + T = (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2)) / (help1) entrop_rho = SVector{ncomponents(equations), real(equations)}(-1.0 * (cv[i] * log(T) - @@ -470,8 +471,8 @@ end rho_v3 = rho * v3 gamma = totalgamma(prim, equations) - rho_e = p / (gamma - 1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + - 0.5 * (B1^2 + B2^2 + B3^2) + rho_e = p / (gamma - 1) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + + 0.5f0 * (B1^2 + B2^2 + B3^2) cons_other = SVector{7, real(equations)}(rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3) @@ -482,9 +483,9 @@ end rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u rho = density(u, equations) gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2)) + 0.5f0 * (B1^2 + B2^2 + B3^2)) return rho * p end @@ -498,7 +499,7 @@ end v3 = rho_v3 / rho v_mag = sqrt(v1^2 + v2^2 + v3^2) gamma = totalgamma(cons, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * v_mag^2 - 0.5 * (B1^2 + B2^2 + B3^2)) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * v_mag^2 - 0.5f0 * (B1^2 + B2^2 + B3^2)) a_square = gamma * p / rho sqrt_rho = sqrt(rho) b1 = B1 / sqrt_rho @@ -506,8 +507,8 @@ end b3 = B3 / sqrt_rho b_square = b1^2 + b2^2 + b3^2 - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) return c_f end diff --git a/src/equations/ideal_glm_mhd_multicomponent_2d.jl b/src/equations/ideal_glm_mhd_multicomponent_2d.jl index a3a50c0485f..51e1dca708e 100644 --- a/src/equations/ideal_glm_mhd_multicomponent_2d.jl +++ b/src/equations/ideal_glm_mhd_multicomponent_2d.jl @@ -129,7 +129,7 @@ function initial_condition_weak_blast_wave(x, t, phi = atan(y_norm, x_norm) sin_phi, cos_phi = sincos(phi) - prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5 ? + prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5f0 ? 2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * @@ -140,9 +140,9 @@ function initial_condition_weak_blast_wave(x, t, 1.1691 for i in eachcomponent(equations)) - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi - p = r > 0.5 ? 1.0 : 1.245 + v1 = r > 0.5f0 ? 0.0 : 0.1882 * cos_phi + v2 = r > 0.5f0 ? 0.0 : 0.1882 * sin_phi + p = r > 0.5f0 ? 1.0 : 1.245 prim_other = SVector{8, real(equations)}(v1, v2, 0.0, p, 1.0, 1.0, 1.0, 0.0) @@ -160,10 +160,10 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * rho * (v1^2 + v2^2 + v3^2) - mag_en = 0.5 * (B1^2 + B2^2 + B3^2) + kin_en = 0.5f0 * rho * (v1^2 + v2^2 + v3^2) + mag_en = 0.5f0 * (B1^2 + B2^2 + B3^2) gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - kin_en - mag_en - 0.5 * psi^2) + p = (gamma - 1) * (rho_e - kin_en - mag_en - 0.5f0 * psi^2) if orientation == 1 f_rho = densities(u, v1, equations) @@ -277,7 +277,7 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 8], u_rr[i + 8]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 8] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 8] + u_rr[i + 8]) for i in eachcomponent(equations)) @@ -287,13 +287,13 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr v3_rr = rho_v3_rr / rho_rr - v1_sq = 0.5 * (v1_ll^2 + v1_rr^2) - v2_sq = 0.5 * (v2_ll^2 + v2_rr^2) - v3_sq = 0.5 * (v3_ll^2 + v3_rr^2) + v1_sq = 0.5f0 * (v1_ll^2 + v1_rr^2) + v2_sq = 0.5f0 * (v2_ll^2 + v2_rr^2) + v3_sq = 0.5f0 * (v3_ll^2 + v3_rr^2) v_sq = v1_sq + v2_sq + v3_sq - B1_sq = 0.5 * (B1_ll^2 + B1_rr^2) - B2_sq = 0.5 * (B2_ll^2 + B2_rr^2) - B3_sq = 0.5 * (B3_ll^2 + B3_rr^2) + B1_sq = 0.5f0 * (B1_ll^2 + B1_rr^2) + B2_sq = 0.5f0 * (B2_ll^2 + B2_rr^2) + B3_sq = 0.5f0 * (B3_ll^2 + B3_rr^2) B_sq = B1_sq + B2_sq + B3_sq vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2 vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 @@ -304,17 +304,17 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, vel_dot_mag_rr = v1_rr * B1_rr + v2_rr * B2_rr + v3_rr * B3_rr # Compute the necessary mean values needed for either direction - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) v_sum = v1_avg + v2_avg + v3_avg - B1_avg = 0.5 * (B1_ll + B1_rr) - B2_avg = 0.5 * (B2_ll + B2_rr) - B3_avg = 0.5 * (B3_ll + B3_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - vel_norm_avg = 0.5 * (vel_norm_ll + vel_norm_rr) - mag_norm_avg = 0.5 * (mag_norm_ll + mag_norm_rr) - vel_dot_mag_avg = 0.5 * (vel_dot_mag_ll + vel_dot_mag_rr) + B1_avg = 0.5f0 * (B1_ll + B1_rr) + B2_avg = 0.5f0 * (B2_ll + B2_rr) + B3_avg = 0.5f0 * (B3_ll + B3_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr) + mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr) + vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr) enth = zero(v_sum) help1_ll = zero(v1_ll) @@ -326,11 +326,11 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, help1_rr += u_rr[i + 8] * cv[i] end - T_ll = (rho_e_ll - 0.5 * rho_ll * (vel_norm_ll) - 0.5 * mag_norm_ll - - 0.5 * psi_ll^2) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (vel_norm_rr) - 0.5 * mag_norm_rr - - 0.5 * psi_rr^2) / help1_rr - T = 0.5 * (1.0 / T_ll + 1.0 / T_rr) + T_ll = (rho_e_ll - 0.5f0 * rho_ll * (vel_norm_ll) - 0.5f0 * mag_norm_ll - + 0.5f0 * psi_ll^2) / help1_ll + T_rr = (rho_e_rr - 0.5f0 * rho_rr * (vel_norm_rr) - 0.5f0 * mag_norm_rr - + 0.5f0 * psi_rr^2) / help1_rr + T = 0.5f0 * (1.0 / T_ll + 1.0 / T_rr) T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) # Calculate fluxes depending on orientation with specific direction averages @@ -343,7 +343,7 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, help1 += f_rho[i] * cv[i] help2 += f_rho[i] end - f1 = help2 * v1_avg + enth / T + 0.5 * mag_norm_avg - B1_avg * B1_avg + f1 = help2 * v1_avg + enth / T + 0.5f0 * mag_norm_avg - B1_avg * B1_avg f2 = help2 * v2_avg - B1_avg * B2_avg f3 = help2 * v3_avg - B1_avg * B3_avg f5 = c_h * psi_avg @@ -351,12 +351,13 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, f7 = v1_avg * B3_avg - v3_avg * B1_avg f8 = c_h * B1_avg # total energy flux is complicated and involves the previous eight components - psi_B1_avg = 0.5 * (B1_ll * psi_ll + B1_rr * psi_rr) - v1_mag_avg = 0.5 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) + psi_B1_avg = 0.5f0 * (B1_ll * psi_ll + B1_rr * psi_rr) + v1_mag_avg = 0.5f0 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) - f4 = (help1 / T_log) - 0.5 * (vel_norm_avg) * (help2) + f1 * v1_avg + + f4 = (help1 / T_log) - 0.5f0 * (vel_norm_avg) * (help2) + f1 * v1_avg + f2 * v2_avg + f3 * v3_avg + - f5 * B1_avg + f6 * B2_avg + f7 * B3_avg + f8 * psi_avg - 0.5 * v1_mag_avg + + f5 * B1_avg + f6 * B2_avg + f7 * B3_avg + f8 * psi_avg - + 0.5f0 * v1_mag_avg + B1_avg * vel_dot_mag_avg - c_h * psi_B1_avg else @@ -367,7 +368,7 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, help2 += f_rho[i] end f1 = help2 * v1_avg - B1_avg * B2_avg - f2 = help2 * v2_avg + enth / T + 0.5 * mag_norm_avg - B2_avg * B2_avg + f2 = help2 * v2_avg + enth / T + 0.5f0 * mag_norm_avg - B2_avg * B2_avg f3 = help2 * v3_avg - B2_avg * B3_avg f5 = v2_avg * B1_avg - v1_avg * B2_avg f6 = c_h * psi_avg @@ -375,12 +376,13 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, f8 = c_h * B2_avg # total energy flux is complicated and involves the previous eight components - psi_B2_avg = 0.5 * (B2_ll * psi_ll + B2_rr * psi_rr) - v2_mag_avg = 0.5 * (v2_ll * mag_norm_ll + v2_rr * mag_norm_rr) + psi_B2_avg = 0.5f0 * (B2_ll * psi_ll + B2_rr * psi_rr) + v2_mag_avg = 0.5f0 * (v2_ll * mag_norm_ll + v2_rr * mag_norm_rr) - f4 = (help1 / T_log) - 0.5 * (vel_norm_avg) * (help2) + f1 * v1_avg + + f4 = (help1 / T_log) - 0.5f0 * (vel_norm_avg) * (help2) + f1 * v1_avg + f2 * v2_avg + f3 * v3_avg + - f5 * B1_avg + f6 * B2_avg + f7 * B3_avg + f8 * psi_avg - 0.5 * v2_mag_avg + + f5 * B1_avg + f6 * B2_avg + f7 * B3_avg + f8 * psi_avg - + 0.5f0 * v2_mag_avg + B2_avg * vel_dot_mag_avg - c_h * psi_B2_avg end @@ -425,20 +427,20 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations. # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_avg = 0.5 * (p_ll + p_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) - magnetic_square_avg = 0.5 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) - inv_gamma_minus_one = 1 / (totalgamma(0.5 * (u_ll + u_rr), equations) - 1) + inv_gamma_minus_one = 1 / (totalgamma(0.5f0 * (u_ll + u_rr), equations) - 1) rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 8], u_rr[i + 8]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 8] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 8] + u_rr[i + 8]) for i in eachcomponent(equations)) @@ -452,18 +454,18 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations. # Calculate fluxes depending on orientation with specific direction averages f2 = f1 * v1_avg + p_avg + magnetic_square_avg - - 0.5 * (B1_ll * B1_rr + B1_rr * B1_ll) - f3 = f1 * v2_avg - 0.5 * (B1_ll * B2_rr + B1_rr * B2_ll) - f4 = f1 * v3_avg - 0.5 * (B1_ll * B3_rr + B1_rr * B3_ll) + 0.5f0 * (B1_ll * B1_rr + B1_rr * B1_ll) + f3 = f1 * v2_avg - 0.5f0 * (B1_ll * B2_rr + B1_rr * B2_ll) + f4 = f1 * v3_avg - 0.5f0 * (B1_ll * B3_rr + B1_rr * B3_ll) #f5 below f6 = f6 = equations.c_h * psi_avg - f7 = 0.5 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) - f8 = 0.5 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) - f9 = equations.c_h * 0.5 * (B1_ll + B1_rr) + f7 = 0.5f0 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) + f8 = 0.5f0 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) + f9 = equations.c_h * 0.5f0 * (B1_ll + B1_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (+p_ll * v1_rr + p_rr * v1_ll + 0.5f0 * (+p_ll * v1_rr + p_rr * v1_ll + (v1_ll * B2_ll * B2_rr + v1_rr * B2_rr * B2_ll) + (v1_ll * B3_ll * B3_rr + v1_rr * B3_rr * B3_ll) - @@ -481,19 +483,19 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations. end # Calculate fluxes depending on orientation with specific direction averages - f2 = f1 * v1_avg - 0.5 * (B2_ll * B1_rr + B2_rr * B1_ll) + f2 = f1 * v1_avg - 0.5f0 * (B2_ll * B1_rr + B2_rr * B1_ll) f3 = f1 * v2_avg + p_avg + magnetic_square_avg - - 0.5 * (B2_ll * B2_rr + B2_rr * B2_ll) - f4 = f1 * v3_avg - 0.5 * (B2_ll * B3_rr + B2_rr * B3_ll) + 0.5f0 * (B2_ll * B2_rr + B2_rr * B2_ll) + f4 = f1 * v3_avg - 0.5f0 * (B2_ll * B3_rr + B2_rr * B3_ll) #f5 below - f6 = 0.5 * (v2_ll * B1_ll - v1_ll * B2_ll + v2_rr * B1_rr - v1_rr * B2_rr) + f6 = 0.5f0 * (v2_ll * B1_ll - v1_ll * B2_ll + v2_rr * B1_rr - v1_rr * B2_rr) f7 = equations.c_h * psi_avg - f8 = 0.5 * (v2_ll * B3_ll - v3_ll * B2_ll + v2_rr * B3_rr - v3_rr * B2_rr) - f9 = equations.c_h * 0.5 * (B2_ll + B2_rr) + f8 = 0.5f0 * (v2_ll * B3_ll - v3_ll * B2_ll + v2_rr * B3_rr - v3_rr * B2_rr) + f9 = equations.c_h * 0.5f0 * (B2_ll + B2_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (+p_ll * v2_rr + p_rr * v2_ll + 0.5f0 * (+p_ll * v2_rr + p_rr * v2_ll + (v2_ll * B1_ll * B1_rr + v2_rr * B1_rr * B1_ll) + (v2_ll * B3_ll * B3_rr + v2_rr * B3_rr * B3_ll) - @@ -550,11 +552,11 @@ end rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u rho = density(u, equations) gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5f0 * (B1^2 + B2^2 + B3^2) - - 0.5 * psi^2) + 0.5f0 * psi^2) return rho * p end @@ -573,8 +575,8 @@ function cons2prim(u, equations::IdealGlmMhdMulticomponentEquations2D) gamma = totalgamma(u, equations) p = (gamma - 1) * - (rho_e - 0.5 * rho * (v1^2 + v2^2 + v3^2) - 0.5 * (B1^2 + B2^2 + B3^2) - - 0.5 * psi^2) + (rho_e - 0.5f0 * rho * (v1^2 + v2^2 + v3^2) - 0.5f0 * (B1^2 + B2^2 + B3^2) - + 0.5f0 * psi^2) prim_other = SVector{8, real(equations)}(v1, v2, v3, p, B1, B2, B3, psi) return vcat(prim_other, prim_rho) end @@ -592,7 +594,7 @@ end v_square = v1^2 + v2^2 + v3^2 gamma = totalgamma(u, equations) p = (gamma - 1) * - (rho_e - 0.5 * rho * v_square - 0.5 * (B1^2 + B2^2 + B3^2) - 0.5 * psi^2) + (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2) - 0.5f0 * psi^2) s = log(p) - gamma * log(rho) rho_p = rho / p @@ -603,7 +605,7 @@ end help1 += u[i + 8] * cv[i] end - T = (rho_e - 0.5 * rho * v_square - 0.5 * (B1^2 + B2^2 + B3^2) - 0.5 * psi^2) / + T = (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2) - 0.5f0 * psi^2) / (help1) entrop_rho = SVector{ncomponents(equations), real(equations)}(-1.0 * @@ -642,8 +644,8 @@ end rho_v3 = rho * v3 gamma = totalgamma(prim, equations) - rho_e = p / (gamma - 1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5 * psi^2 + rho_e = p / (gamma - 1) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + + 0.5f0 * (B1^2 + B2^2 + B3^2) + 0.5f0 * psi^2 cons_other = SVector{8, real(equations)}(rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi) @@ -662,7 +664,7 @@ end v_mag = sqrt(v1^2 + v2^2 + v3^2) gamma = totalgamma(cons, equations) p = (gamma - 1) * - (rho_e - 0.5 * rho * v_mag^2 - 0.5 * (B1^2 + B2^2 + B3^2) - 0.5 * psi^2) + (rho_e - 0.5f0 * rho * v_mag^2 - 0.5f0 * (B1^2 + B2^2 + B3^2) - 0.5f0 * psi^2) a_square = gamma * p / rho sqrt_rho = sqrt(rho) b1 = B1 / sqrt_rho @@ -670,11 +672,11 @@ end b3 = B3 / sqrt_rho b_square = b1^2 + b2^2 + b3^2 if direction == 1 # x-direction - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) else - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b2^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b2^2)) end return c_f end diff --git a/src/equations/inviscid_burgers_1d.jl b/src/equations/inviscid_burgers_1d.jl index 130196a4929..b9b0e67d6af 100644 --- a/src/equations/inviscid_burgers_1d.jl +++ b/src/equations/inviscid_burgers_1d.jl @@ -69,7 +69,7 @@ end # Calculate 1D flux in for a single point @inline function flux(u, orientation::Integer, equation::InviscidBurgersEquation1D) - return SVector(0.5 * u[1]^2) + return SVector(0.5f0 * u[1]^2) end # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation @@ -111,7 +111,7 @@ function flux_godunov(u_ll, u_rr, orientation, equation::InviscidBurgersEquation u_L = u_ll[1] u_R = u_rr[1] - return SVector(0.5 * max(max(u_L, zero(u_L))^2, min(u_R, zero(u_R))^2)) + return SVector(0.5f0 * max(max(u_L, zero(u_L))^2, min(u_R, zero(u_R))^2)) end # See https://metaphor.ethz.ch/x/2019/hs/401-4671-00L/literature/mishra_hyperbolic_pdes.pdf , @@ -121,7 +121,7 @@ function flux_engquist_osher(u_ll, u_rr, orientation, u_L = u_ll[1] u_R = u_rr[1] - return SVector(0.5 * (max(u_L, zero(u_L))^2 + min(u_R, zero(u_R))^2)) + return SVector(0.5f0 * (max(u_L, zero(u_L))^2 + min(u_R, zero(u_R))^2)) end """ @@ -151,16 +151,16 @@ end @inline function splitting_lax_friedrichs(u, ::Val{:plus}, orientation::Integer, equations::InviscidBurgersEquation1D) - f = 0.5 * u[1]^2 + f = 0.5f0 * u[1]^2 lambda = abs(u[1]) - return SVector(0.5 * (f + lambda * u[1])) + return SVector(0.5f0 * (f + lambda * u[1])) end @inline function splitting_lax_friedrichs(u, ::Val{:minus}, orientation::Integer, equations::InviscidBurgersEquation1D) - f = 0.5 * u[1]^2 + f = 0.5f0 * u[1]^2 lambda = abs(u[1]) - return SVector(0.5 * (f - lambda * u[1])) + return SVector(0.5f0 * (f - lambda * u[1])) end # Convert conservative variables to primitive @@ -171,11 +171,11 @@ end @inline entropy2cons(u, equation::InviscidBurgersEquation1D) = u # Calculate entropy for a conservative state `cons` -@inline entropy(u::Real, ::InviscidBurgersEquation1D) = 0.5 * u^2 +@inline entropy(u::Real, ::InviscidBurgersEquation1D) = 0.5f0 * u^2 @inline entropy(u, equation::InviscidBurgersEquation1D) = entropy(u[1], equation) # Calculate total energy for a conservative state `cons` -@inline energy_total(u::Real, ::InviscidBurgersEquation1D) = 0.5 * u^2 +@inline energy_total(u::Real, ::InviscidBurgersEquation1D) = 0.5f0 * u^2 @inline function energy_total(u, equation::InviscidBurgersEquation1D) energy_total(u[1], equation) end From 772c63f3d50b111087bfac6e2f28c6f0fd5101f4 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 1 Jul 2024 07:29:49 -1000 Subject: [PATCH 2/6] continue --- src/equations/ideal_glm_mhd_multicomponent_1d.jl | 10 ++++++---- src/equations/ideal_glm_mhd_multicomponent_2d.jl | 6 +++--- src/equations/inviscid_burgers_1d.jl | 4 ++-- 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/src/equations/ideal_glm_mhd_multicomponent_1d.jl b/src/equations/ideal_glm_mhd_multicomponent_1d.jl index 7ed3fef9d82..18ae51e6daf 100644 --- a/src/equations/ideal_glm_mhd_multicomponent_1d.jl +++ b/src/equations/ideal_glm_mhd_multicomponent_1d.jl @@ -83,18 +83,20 @@ function initial_condition_convergence_test(x, t, # smooth Alfvén wave test from Derigs et al. FLASH (2016) # domain must be set to [0, 1], γ = 5/3 - rho = 1.0 + RealT = eltype(x) + rho = 1 prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * rho for i in eachcomponent(equations)) - v1 = 0.0 - si, co = sincos(2 * pi * x[1]) + v1 = 0 + # TODO: sincospi + si, co = sincos(2 * convert(RealT, pi) * x[1]) v2 = 0.1 * si v3 = 0.1 * co p = 0.1 - B1 = 1.0 + B1 = 1 B2 = v2 B3 = v3 prim_other = SVector{7, real(equations)}(v1, v2, v3, p, B1, B2, B3) diff --git a/src/equations/ideal_glm_mhd_multicomponent_2d.jl b/src/equations/ideal_glm_mhd_multicomponent_2d.jl index 51e1dca708e..2ac62235fea 100644 --- a/src/equations/ideal_glm_mhd_multicomponent_2d.jl +++ b/src/equations/ideal_glm_mhd_multicomponent_2d.jl @@ -88,9 +88,9 @@ function initial_condition_convergence_test(x, t, equations::IdealGlmMhdMulticomponentEquations2D) # smooth Alfvén wave test from Derigs et al. FLASH (2016) # domain must be set to [0, 1/cos(α)] x [0, 1/sin(α)], γ = 5/3 - alpha = 0.25 * pi + alpha = 0.25 * convert(RealT, pi) x_perp = x[1] * cos(alpha) + x[2] * sin(alpha) - B_perp = 0.1 * sin(2.0 * pi * x_perp) + B_perp = 0.1 * sinpi(2 * x_perp) rho = 1 prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / (1 - @@ -99,7 +99,7 @@ function initial_condition_convergence_test(x, t, for i in eachcomponent(equations)) v1 = -B_perp * sin(alpha) v2 = B_perp * cos(alpha) - v3 = 0.1 * cos(2.0 * pi * x_perp) + v3 = 0.1 * cospi(2 * x_perp) p = 0.1 B1 = cos(alpha) + v1 B2 = sin(alpha) + v2 diff --git a/src/equations/inviscid_burgers_1d.jl b/src/equations/inviscid_burgers_1d.jl index b9b0e67d6af..4091663d729 100644 --- a/src/equations/inviscid_burgers_1d.jl +++ b/src/equations/inviscid_burgers_1d.jl @@ -39,7 +39,7 @@ function initial_condition_convergence_test(x, t, equation::InviscidBurgersEquat A = 1.0 L = 1 f = 1 / L - omega = 2 * pi * f + omega = 2 * convert(RealT, pi) * f scalar = c + A * sin(omega * (x[1] - t)) return SVector(scalar) @@ -58,7 +58,7 @@ Source terms used for convergence tests in combination with A = 1.0 L = 1 f = 1 / L - omega = 2 * pi * f + omega = 2 * convert(RealT, pi) * f du = omega * A * cos(omega * (x[1] - t)) * (c - 1 + A * sin(omega * (x[1] - t))) return SVector(du) From 7e148e9f38cfaab3df371f7bc6565ce1783db024 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 1 Jul 2024 21:14:31 -1000 Subject: [PATCH 3/6] complete equations --- .../compressible_euler_multicomponent_1d.jl | 2 +- .../compressible_euler_multicomponent_2d.jl | 2 +- .../ideal_glm_mhd_multicomponent_1d.jl | 76 +++++++++-------- .../ideal_glm_mhd_multicomponent_2d.jl | 83 ++++++++++--------- src/equations/inviscid_burgers_1d.jl | 17 ++-- 5 files changed, 100 insertions(+), 80 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl index d4f6b421f7c..da434579f76 100644 --- a/src/equations/compressible_euler_multicomponent_1d.jl +++ b/src/equations/compressible_euler_multicomponent_1d.jl @@ -206,7 +206,7 @@ function initial_condition_weak_blast_wave(x, t, 2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * - 1 : + one(RealT) : 2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 3473f887336..2424ad0301c 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -224,7 +224,7 @@ function initial_condition_weak_blast_wave(x, t, 2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * - 1 : + one(RealT) : 2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * diff --git a/src/equations/ideal_glm_mhd_multicomponent_1d.jl b/src/equations/ideal_glm_mhd_multicomponent_1d.jl index 18ae51e6daf..b2ed06e53ea 100644 --- a/src/equations/ideal_glm_mhd_multicomponent_1d.jl +++ b/src/equations/ideal_glm_mhd_multicomponent_1d.jl @@ -93,13 +93,14 @@ function initial_condition_convergence_test(x, t, v1 = 0 # TODO: sincospi si, co = sincos(2 * convert(RealT, pi) * x[1]) - v2 = 0.1 * si - v3 = 0.1 * co - p = 0.1 + v2 = convert(RealT, 0.1) * si + v3 = convert(RealT, 0.1) * co + p = convert(RealT, 0.1) B1 = 1 B2 = v2 B3 = v3 - prim_other = SVector{7, real(equations)}(v1, v2, v3, p, B1, B2, B3) + prim_other = SVector(v1, v2, v3, p, B1, B2, B3) + return prim2cons(vcat(prim_other, prim_rho), equations) end @@ -116,6 +117,7 @@ function initial_condition_weak_blast_wave(x, t, # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # Same discontinuity in the velocities but with magnetic fields # Set up polar coordinates + RealT = eltype(x) inicenter = (0) x_norm = x[1] - inicenter[1] r = sqrt(x_norm^2) @@ -123,24 +125,24 @@ function initial_condition_weak_blast_wave(x, t, # Calculate primitive variables if r > 0.5f0 - rho = 1.0 + rho = one(RealT) prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * rho for i in eachcomponent(equations)) else - rho = 1.1691 + rho = convert(RealT, 1.1691) prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * rho for i in eachcomponent(equations)) end - v1 = r > 0.5f0 ? 0.0 : 0.1882 * cos(phi) - p = r > 0.5f0 ? 1.0 : 1.245 + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi) + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) - prim_other = SVector{7, real(equations)}(v1, 0.0, 0.0, p, 1.0, 1.0, 1.0) + prim_other = SVector(v1, 0, 0, p, 1, 1, 1) return prim2cons(vcat(prim_other, prim_rho), equations) end @@ -166,11 +168,11 @@ end f3 = rho_v1 * v3 - B1 * B3 f4 = (kin_en + gamma * p / (gamma - 1) + 2 * mag_en) * v1 - B1 * (v1 * B1 + v2 * B2 + v3 * B3) - f5 = 0.0 + f5 = 0 f6 = v1 * B2 - v2 * B1 f7 = v1 * B3 - v3 * B1 - f_other = SVector{7, real(equations)}(f1, f2, f3, f4, f5, f6, f7) + f_other = SVector(f1, f2, f3, f4, f5, f6, f7) return vcat(f_other, f_rho) end @@ -230,9 +232,10 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr) vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr) - enth = zero(v_sum) - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) + RealT = eltype(u_ll) + enth = zero(RealT) + help1_ll = zero(RealT) + help1_rr = zero(RealT) for i in eachcomponent(equations) enth += rhok_avg[i] * gas_constants[i] @@ -242,12 +245,12 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, T_ll = (rho_e_ll - 0.5f0 * rho_ll * (vel_norm_ll) - 0.5f0 * mag_norm_ll) / help1_ll T_rr = (rho_e_rr - 0.5f0 * rho_rr * (vel_norm_rr) - 0.5f0 * mag_norm_rr) / help1_rr - T = 0.5f0 * (1.0 / T_ll + 1.0 / T_rr) - T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) + T = 0.5f0 * (1 / T_ll + 1 / T_rr) + T_log = ln_mean(1 / T_ll, 1 / T_rr) # Calculate fluxes depending on orientation with specific direction averages - help1 = zero(T_ll) - help2 = zero(T_rr) + help1 = zero(RealT) + help2 = zero(RealT) f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg for i in eachcomponent(equations)) @@ -258,7 +261,7 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, f1 = help2 * v1_avg + enth / T + 0.5f0 * mag_norm_avg - B1_avg * B1_avg f2 = help2 * v2_avg - B1_avg * B2_avg f3 = help2 * v3_avg - B1_avg * B3_avg - f5 = 0.0 + f5 = 0 f6 = v1_avg * B2_avg - v2_avg * B1_avg f7 = v1_avg * B3_avg - v3_avg * B1_avg @@ -271,7 +274,7 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, f5 * B1_avg + f6 * B2_avg + f7 * B3_avg - 0.5f0 * v1_mag_avg + B1_avg * vel_dot_mag_avg - f_other = SVector{7, real(equations)}(f1, f2, f3, f4, f5, f6, f7) + f_other = SVector(f1, f2, f3, f4, f5, f6, f7) return vcat(f_other, f_rho) end @@ -328,7 +331,8 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations. u_rr[i + 7]) for i in eachcomponent(equations)) - f1 = zero(rho_ll) + RealT = eltype(u_ll) + f1 = zero(RealT) f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg for i in eachcomponent(equations)) for i in eachcomponent(equations) @@ -340,8 +344,8 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations. 0.5f0 * (B1_ll * B1_rr + B1_rr * B1_ll) f3 = f1 * v2_avg - 0.5f0 * (B1_ll * B2_rr + B1_rr * B2_ll) f4 = f1 * v3_avg - 0.5f0 * (B1_ll * B3_rr + B1_rr * B3_ll) - #f5 below - f6 = 0.0 + # f5 below + f6 = 0 f7 = 0.5f0 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) f8 = 0.5f0 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) # total energy flux is complicated and involves the previous components @@ -355,7 +359,7 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations. - (v3_ll * B1_ll * B3_rr + v3_rr * B1_rr * B3_ll))) - f_other = SVector{7, real(equations)}(f2, f3, f4, f5, f6, f7, f8) + f_other = SVector(f2, f3, f4, f5, f6, f7, f8) return vcat(f_other, f_rho) end @@ -409,7 +413,8 @@ function cons2prim(u, equations::IdealGlmMhdMulticomponentEquations1D) p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2 + v2^2 + v3^2) - 0.5f0 * (B1^2 + B2^2 + B3^2)) - prim_other = SVector{7, real(equations)}(v1, v2, v3, p, B1, B2, B3) + prim_other = SVector(v1, v2, v3, p, B1, B2, B3) + return vcat(prim_other, prim_rho) end @@ -430,7 +435,8 @@ end rho_p = rho / p # Multicomponent stuff - help1 = zero(v1) + RealT = eltype(u) + help1 = zero(RealT) for i in eachcomponent(equations) help1 += u[i + 7] * cv[i] @@ -438,7 +444,7 @@ end T = (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2)) / (help1) - entrop_rho = SVector{ncomponents(equations), real(equations)}(-1.0 * + entrop_rho = SVector{ncomponents(equations), real(equations)}(-1 * (cv[i] * log(T) - gas_constants[i] * log(u[i + 7])) + @@ -450,12 +456,12 @@ end w1 = v1 / T w2 = v2 / T w3 = v3 / T - w4 = -1.0 / T + w4 = -1 / T w5 = B1 / T w6 = B2 / T w7 = B3 / T - entrop_other = SVector{7, real(equations)}(w1, w2, w3, w4, w5, w6, w7) + entrop_other = SVector(w1, w2, w3, w4, w5, w6, w7) return vcat(entrop_other, entrop_rho) end @@ -476,7 +482,7 @@ end rho_e = p / (gamma - 1) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + 0.5f0 * (B1^2 + B2^2 + B3^2) - cons_other = SVector{7, real(equations)}(rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3) + cons_other = SVector(rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3) return vcat(cons_other, cons_rho) end @@ -510,13 +516,14 @@ end b_square = b1^2 + b2^2 + b3^2 c_f = sqrt(0.5f0 * (a_square + b_square) + - 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b1^2)) return c_f end @inline function density(u, equations::IdealGlmMhdMulticomponentEquations1D) - rho = zero(u[1]) + RealT = eltype(u) + rho = zero(RealT) for i in eachcomponent(equations) rho += u[i + 7] @@ -528,8 +535,9 @@ end @inline function totalgamma(u, equations::IdealGlmMhdMulticomponentEquations1D) @unpack cv, gammas = equations - help1 = zero(u[1]) - help2 = zero(u[1]) + RealT = eltype(u) + help1 = zero(RealT) + help2 = zero(RealT) for i in eachcomponent(equations) help1 += u[i + 7] * cv[i] * gammas[i] diff --git a/src/equations/ideal_glm_mhd_multicomponent_2d.jl b/src/equations/ideal_glm_mhd_multicomponent_2d.jl index 2ac62235fea..3aab048bd99 100644 --- a/src/equations/ideal_glm_mhd_multicomponent_2d.jl +++ b/src/equations/ideal_glm_mhd_multicomponent_2d.jl @@ -88,24 +88,27 @@ function initial_condition_convergence_test(x, t, equations::IdealGlmMhdMulticomponentEquations2D) # smooth Alfvén wave test from Derigs et al. FLASH (2016) # domain must be set to [0, 1/cos(α)] x [0, 1/sin(α)], γ = 5/3 - alpha = 0.25 * convert(RealT, pi) + RealT = eltype(x) + alpha = 0.25f0 * convert(RealT, pi) x_perp = x[1] * cos(alpha) + x[2] * sin(alpha) - B_perp = 0.1 * sinpi(2 * x_perp) + B_perp = convert(RealT, 0.1) * sinpi(2 * x_perp) rho = 1 prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * rho for i in eachcomponent(equations)) + v1 = -B_perp * sin(alpha) v2 = B_perp * cos(alpha) - v3 = 0.1 * cospi(2 * x_perp) - p = 0.1 + v3 = convert(RealT, 0.1) * cospi(2 * x_perp) + p = convert(RealT, 0.1) B1 = cos(alpha) + v1 B2 = sin(alpha) + v2 B3 = v3 - psi = 0.0 - prim_other = SVector{8, real(equations)}(v1, v2, v3, p, B1, B2, B3, psi) + psi = 0 + prim_other = SVector(v1, v2, v3, p, B1, B2, B3, psi) + return prim2cons(vcat(prim_other, prim_rho), equations) end @@ -122,7 +125,8 @@ function initial_condition_weak_blast_wave(x, t, # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # Same discontinuity in the velocities but with magnetic fields # Set up polar coordinates - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) @@ -133,18 +137,18 @@ function initial_condition_weak_blast_wave(x, t, 2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * - 1.0 : + one(RealT) : 2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * - 1.1691 + convert(RealT, 1.1691) for i in eachcomponent(equations)) - v1 = r > 0.5f0 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5f0 ? 0.0 : 0.1882 * sin_phi - p = r > 0.5f0 ? 1.0 : 1.245 + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) - prim_other = SVector{8, real(equations)}(v1, v2, 0.0, p, 1.0, 1.0, 1.0, 0.0) + prim_other = SVector(v1, v2, 0, p, 1, 1, 1, 0) return prim2cons(vcat(prim_other, prim_rho), equations) end @@ -189,7 +193,7 @@ end f8 = c_h * B2 end - f_other = SVector{8, real(equations)}(f1, f2, f3, f4, f5, f6, f7, f8) + f_other = SVector(f1, f2, f3, f4, f5, f6, f7, f8) return vcat(f_other, f_rho) end @@ -316,9 +320,10 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr) vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr) - enth = zero(v_sum) - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) + RealT = eltype(u_ll) + enth = zero(RealT) + help1_ll = zero(RealT) + help1_rr = zero(RealT) for i in eachcomponent(equations) enth += rhok_avg[i] * gas_constants[i] @@ -330,12 +335,12 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, 0.5f0 * psi_ll^2) / help1_ll T_rr = (rho_e_rr - 0.5f0 * rho_rr * (vel_norm_rr) - 0.5f0 * mag_norm_rr - 0.5f0 * psi_rr^2) / help1_rr - T = 0.5f0 * (1.0 / T_ll + 1.0 / T_rr) - T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) + T = 0.5f0 * (1 / T_ll + 1 / T_rr) + T_log = ln_mean(1 / T_ll, 1 / T_rr) # Calculate fluxes depending on orientation with specific direction averages - help1 = zero(T_ll) - help2 = zero(T_rr) + help1 = zero(RealT) + help2 = zero(RealT) if orientation == 1 f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg for i in eachcomponent(equations)) @@ -386,7 +391,7 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, B2_avg * vel_dot_mag_avg - c_h * psi_B2_avg end - f_other = SVector{8, real(equations)}(f1, f2, f3, f4, f5, f6, f7, f8) + f_other = SVector(f1, f2, f3, f4, f5, f6, f7, f8) return vcat(f_other, f_rho) end @@ -444,8 +449,9 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations. u_rr[i + 8]) for i in eachcomponent(equations)) + RealT = eltype(u_ll) if orientation == 1 - f1 = zero(rho_ll) + f1 = zero(RealT) f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg for i in eachcomponent(equations)) for i in eachcomponent(equations) @@ -457,7 +463,7 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations. 0.5f0 * (B1_ll * B1_rr + B1_rr * B1_ll) f3 = f1 * v2_avg - 0.5f0 * (B1_ll * B2_rr + B1_rr * B2_ll) f4 = f1 * v3_avg - 0.5f0 * (B1_ll * B3_rr + B1_rr * B3_ll) - #f5 below + # f5 below f6 = f6 = equations.c_h * psi_avg f7 = 0.5f0 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) f8 = 0.5f0 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) @@ -475,7 +481,7 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations. + equations.c_h * (B1_ll * psi_rr + B1_rr * psi_ll))) else - f1 = zero(rho_ll) + f1 = zero(RealT) f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v2_avg for i in eachcomponent(equations)) for i in eachcomponent(equations) @@ -506,7 +512,7 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations. equations.c_h * (B2_ll * psi_rr + B2_rr * psi_ll))) end - f_other = SVector{8, real(equations)}(f2, f3, f4, f5, f6, f7, f8, f9) + f_other = SVector(f2, f3, f4, f5, f6, f7, f8, f9) return vcat(f_other, f_rho) end @@ -577,7 +583,8 @@ function cons2prim(u, equations::IdealGlmMhdMulticomponentEquations2D) p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2 + v2^2 + v3^2) - 0.5f0 * (B1^2 + B2^2 + B3^2) - 0.5f0 * psi^2) - prim_other = SVector{8, real(equations)}(v1, v2, v3, p, B1, B2, B3, psi) + prim_other = SVector(v1, v2, v3, p, B1, B2, B3, psi) + return vcat(prim_other, prim_rho) end @@ -608,7 +615,7 @@ end T = (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2) - 0.5f0 * psi^2) / (help1) - entrop_rho = SVector{ncomponents(equations), real(equations)}(-1.0 * + entrop_rho = SVector{ncomponents(equations), real(equations)}(-1 * (cv[i] * log(T) - gas_constants[i] * log(u[i + 8])) + @@ -620,13 +627,13 @@ end w1 = v1 / T w2 = v2 / T w3 = v3 / T - w4 = -1.0 / T + w4 = -1 / T w5 = B1 / T w6 = B2 / T w7 = B3 / T w8 = psi / T - entrop_other = SVector{8, real(equations)}(w1, w2, w3, w4, w5, w6, w7, w8) + entrop_other = SVector(w1, w2, w3, w4, w5, w6, w7, w8) return vcat(entrop_other, entrop_rho) end @@ -647,8 +654,8 @@ end rho_e = p / (gamma - 1) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + 0.5f0 * (B1^2 + B2^2 + B3^2) + 0.5f0 * psi^2 - cons_other = SVector{8, real(equations)}(rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, - psi) + cons_other = SVector(rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, + psi) return vcat(cons_other, cons_rho) end @@ -673,16 +680,17 @@ end b_square = b1^2 + b2^2 + b3^2 if direction == 1 # x-direction c_f = sqrt(0.5f0 * (a_square + b_square) + - 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b1^2)) else c_f = sqrt(0.5f0 * (a_square + b_square) + - 0.5f0 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b2^2)) + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b2^2)) end return c_f end @inline function density(u, equations::IdealGlmMhdMulticomponentEquations2D) - rho = zero(u[1]) + RealT = eltype(u) + rho = zero(RealT) for i in eachcomponent(equations) rho += u[i + 8] @@ -694,8 +702,9 @@ end @inline function totalgamma(u, equations::IdealGlmMhdMulticomponentEquations2D) @unpack cv, gammas = equations - help1 = zero(u[1]) - help2 = zero(u[1]) + RealT = eltype(u) + help1 = zero(RealT) + help2 = zero(RealT) for i in eachcomponent(equations) help1 += u[i + 8] * cv[i] * gammas[i] diff --git a/src/equations/inviscid_burgers_1d.jl b/src/equations/inviscid_burgers_1d.jl index 4091663d729..8528f4e56db 100644 --- a/src/equations/inviscid_burgers_1d.jl +++ b/src/equations/inviscid_burgers_1d.jl @@ -26,7 +26,8 @@ varnames(::typeof(cons2prim), ::InviscidBurgersEquation1D) = ("scalar",) A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equation::InviscidBurgersEquation1D) - return SVector(2.0) + RealT = eltype(x) + return SVector(RealT(2)) end """ @@ -35,8 +36,9 @@ end A smooth initial condition used for convergence tests. """ function initial_condition_convergence_test(x, t, equation::InviscidBurgersEquation1D) - c = 2.0 - A = 1.0 + RealT = eltype(x) + c = 2 + A = 1 L = 1 f = 1 / L omega = 2 * convert(RealT, pi) * f @@ -54,8 +56,9 @@ Source terms used for convergence tests in combination with @inline function source_terms_convergence_test(u, x, t, equations::InviscidBurgersEquation1D) # Same settings as in `initial_condition` - c = 2.0 - A = 1.0 + RealT = eltype(x) + c = 2 + A = 1 L = 1 f = 1 / L omega = 2 * convert(RealT, pi) * f @@ -111,7 +114,7 @@ function flux_godunov(u_ll, u_rr, orientation, equation::InviscidBurgersEquation u_L = u_ll[1] u_R = u_rr[1] - return SVector(0.5f0 * max(max(u_L, zero(u_L))^2, min(u_R, zero(u_R))^2)) + return SVector(0.5f0 * max(max(u_L, 0)^2, min(u_R, 0)^2)) end # See https://metaphor.ethz.ch/x/2019/hs/401-4671-00L/literature/mishra_hyperbolic_pdes.pdf , @@ -121,7 +124,7 @@ function flux_engquist_osher(u_ll, u_rr, orientation, u_L = u_ll[1] u_R = u_rr[1] - return SVector(0.5f0 * (max(u_L, zero(u_L))^2 + min(u_R, zero(u_R))^2)) + return SVector(0.5f0 * (max(u_L, 0)^2 + min(u_R, 0)^2)) end """ From 24f446489bd1ecddc5c26fe72eaeef9b90439283 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 4 Jul 2024 13:44:51 -1000 Subject: [PATCH 4/6] unit test 1D --- test/test_type.jl | 58 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/test/test_type.jl b/test/test_type.jl index fc9a9561e6c..4f0d037227e 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1351,6 +1351,64 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Ideal Glm Mhd Multicomponent 1D" begin + for RealT in (Float32, Float64) + gammas = (RealT(2), RealT(2)) + gas_constants = (RealT(2), RealT(2)) + equations = @inferred IdealGlmMhdMulticomponentEquations1D(gammas = gammas, + gas_constants = gas_constants) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT)) + orientation = 1 + directions = [1, 2] + + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, + equations)) == + RealT + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + orientation, + equations)) == RealT + else + @test eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, orientation, + equations)) == RealT + end + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == + RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + @test typeof(@inferred Trixi.totalgamma(u, equations)) == RealT + + for direction in directions + @test typeof(Trixi.calc_fast_wavespeed(cons, direction, equations)) == RealT + end + end + end + @timed_testset "Linear Scalar Advection 1D" begin for RealT in (Float32, Float64) equations = @inferred LinearScalarAdvectionEquation1D(RealT(1)) From 3655e1d076eeb8166d403efdb9f3a62d26563e9d Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 4 Jul 2024 13:59:06 -1000 Subject: [PATCH 5/6] unit test 2D --- test/test_type.jl | 66 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/test/test_type.jl b/test/test_type.jl index 4f0d037227e..a46c2a92bde 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1409,6 +1409,72 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Ideal Glm Mhd Multicomponent 2D" begin + for RealT in (Float32, Float64) + gammas = (RealT(2), RealT(2)) + gas_constants = (RealT(2), RealT(2)) + equations = @inferred IdealGlmMhdMulticomponentEquations2D(gammas = gammas, + gas_constants = gas_constants) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT)) + orientations = [1, 2] + directions = [1, 2, 3, 4] + + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, orientation, + equations)) == + RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, + equations)) == RealT + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + orientation, + equations)) == + RealT + else + @test eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, + equations)) == + RealT + @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, orientation, + equations)) == RealT + end + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == + RealT + end + + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + @test typeof(@inferred Trixi.totalgamma(u, equations)) == RealT + + for direction in directions + @test typeof(Trixi.calc_fast_wavespeed(cons, direction, equations)) == RealT + end + end + end + @timed_testset "Linear Scalar Advection 1D" begin for RealT in (Float32, Float64) equations = @inferred LinearScalarAdvectionEquation1D(RealT(1)) From 7160962013e676a75c99a7e805da324e70d3892f Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 4 Jul 2024 14:26:03 -1000 Subject: [PATCH 6/6] complete unit tests --- src/equations/inviscid_burgers_1d.jl | 4 +-- test/test_type.jl | 42 ++++++++++++++++++++++++++++ 2 files changed, 44 insertions(+), 2 deletions(-) diff --git a/src/equations/inviscid_burgers_1d.jl b/src/equations/inviscid_burgers_1d.jl index 8528f4e56db..afc6f9999c7 100644 --- a/src/equations/inviscid_burgers_1d.jl +++ b/src/equations/inviscid_burgers_1d.jl @@ -40,7 +40,7 @@ function initial_condition_convergence_test(x, t, equation::InviscidBurgersEquat c = 2 A = 1 L = 1 - f = 1 / L + f = 1.0f0 / L omega = 2 * convert(RealT, pi) * f scalar = c + A * sin(omega * (x[1] - t)) @@ -60,7 +60,7 @@ Source terms used for convergence tests in combination with c = 2 A = 1 L = 1 - f = 1 / L + f = 1.0f0 / L omega = 2 * convert(RealT, pi) * f du = omega * A * cos(omega * (x[1] - t)) * (c - 1 + A * sin(omega * (x[1] - t))) diff --git a/test/test_type.jl b/test/test_type.jl index a46c2a92bde..c45a750c07f 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1475,6 +1475,48 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Inviscid Burgers 1D" begin + for RealT in (Float32, Float64) + equations = @inferred InviscidBurgersEquation1D() + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = SVector(one(RealT)) + orientation = 1 + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_ec(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(@inferred flux_godunov(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred Trixi.flux_engquist_osher(u_ll, u_rr, orientation, + equations)) == + RealT + + @test eltype(eltype(@inferred splitting_lax_friedrichs(u, orientation, + equations))) == + RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == + RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test eltype(@inferred entropy2cons(u, equations)) == RealT + @test typeof(@inferred entropy(u, equations)) == RealT + @test typeof(@inferred energy_total(u, equations)) == RealT + end + end + @timed_testset "Linear Scalar Advection 1D" begin for RealT in (Float32, Float64) equations = @inferred LinearScalarAdvectionEquation1D(RealT(1))