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 dad7c27e86c..b2ed06e53ea 100644 --- a/src/equations/ideal_glm_mhd_multicomponent_1d.jl +++ b/src/equations/ideal_glm_mhd_multicomponent_1d.jl @@ -83,21 +83,24 @@ 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]) - v2 = 0.1 * si - v3 = 0.1 * co - p = 0.1 - B1 = 1.0 + v1 = 0 + # TODO: sincospi + si, co = sincos(2 * convert(RealT, pi) * x[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 @@ -114,31 +117,32 @@ 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) phi = atan(x_norm) # Calculate primitive variables - if r > 0.5 - rho = 1.0 + if r > 0.5f0 + 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.5 ? 0.0 : 0.1882 * cos(phi) - p = r > 0.5 ? 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 @@ -153,8 +157,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) @@ -164,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 @@ -198,7 +202,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,20 +221,21 @@ 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) - - enth = zero(v_sum) - help1_ll = zero(v1_ll) - help1_rr = zero(v1_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) + + 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] @@ -238,14 +243,14 @@ 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_log = ln_mean(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 / 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)) @@ -253,22 +258,23 @@ 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 + f5 = 0 f6 = v1_avg * B2_avg - v2_avg * B1_avg 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) + f_other = SVector(f1, f2, f3, f4, f5, f6, f7) return vcat(f_other, f_rho) end @@ -309,23 +315,24 @@ 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)) - 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) @@ -334,17 +341,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) - #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) + 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 + 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) - @@ -352,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 @@ -405,8 +412,9 @@ 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)) - prim_other = SVector{7, real(equations)}(v1, v2, v3, p, B1, B2, B3) + (rho_e - 0.5f0 * rho * (v1^2 + v2^2 + v3^2) - 0.5f0 * (B1^2 + B2^2 + B3^2)) + prim_other = SVector(v1, v2, v3, p, B1, B2, B3) + return vcat(prim_other, prim_rho) end @@ -422,20 +430,21 @@ 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 # Multicomponent stuff - help1 = zero(v1) + RealT = eltype(u) + help1 = zero(RealT) for i in eachcomponent(equations) 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 * + entrop_rho = SVector{ncomponents(equations), real(equations)}(-1 * (cv[i] * log(T) - gas_constants[i] * log(u[i + 7])) + @@ -447,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 @@ -470,10 +479,10 @@ 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) + cons_other = SVector(rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3) return vcat(cons_other, cons_rho) end @@ -482,9 +491,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 +507,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,14 +515,15 @@ 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 * 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] @@ -525,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 a3a50c0485f..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 * pi + RealT = eltype(x) + alpha = 0.25f0 * 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 = 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 * cos(2.0 * pi * 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,29 +125,30 @@ 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) 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)) * - 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.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 ? 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 @@ -160,10 +164,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) @@ -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 @@ -277,7 +281,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 +291,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,21 +308,22 @@ 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) - - enth = zero(v_sum) - help1_ll = zero(v1_ll) - help1_rr = zero(v1_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) + + 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] @@ -326,16 +331,16 @@ 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_log = ln_mean(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 / 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)) @@ -343,7 +348,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 +356,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 +373,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,16 +381,17 @@ 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 - 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 @@ -425,25 +432,26 @@ 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)) + 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) @@ -452,18 +460,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) - #f5 below + 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) - @@ -473,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) @@ -481,19 +489,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) - @@ -504,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 @@ -550,11 +558,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,9 +581,10 @@ 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) - prim_other = SVector{8, real(equations)}(v1, v2, v3, p, B1, B2, B3, psi) + (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(v1, v2, v3, p, B1, B2, B3, psi) + return vcat(prim_other, prim_rho) end @@ -592,7 +601,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,10 +612,10 @@ 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 * + entrop_rho = SVector{ncomponents(equations), real(equations)}(-1 * (cv[i] * log(T) - gas_constants[i] * log(u[i + 8])) + @@ -618,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 @@ -642,11 +651,11 @@ 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) + cons_other = SVector(rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, + psi) return vcat(cons_other, cons_rho) end @@ -662,7 +671,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,17 +679,18 @@ 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 * 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 * 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] @@ -692,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 130196a4929..afc6f9999c7 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,11 +36,12 @@ 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 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f scalar = c + A * sin(omega * (x[1] - t)) return SVector(scalar) @@ -54,11 +56,12 @@ 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 * pi * f + 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))) return SVector(du) @@ -69,7 +72,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 +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.5 * 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.5 * (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 """ @@ -151,16 +154,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 +174,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 diff --git a/test/test_type.jl b/test/test_type.jl index fc9a9561e6c..c45a750c07f 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1351,6 +1351,172 @@ 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 "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 "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))