From 7b57dd52dbee224bdeba352cbda44178cab89e2f Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Fri, 14 Jun 2024 03:57:10 -1000 Subject: [PATCH] Add numerical support of other real types (continue) (#1947) * start main continue * rerun CI tests * test CI tests * complete the first * compressible euler quasi 1d * compressible euler multicomponent 1d * compressible euler multicomponent 2d * complete unit tests * fix comments * complete * apply strict type check Co-authored-by: Hendrik Ranocha * cover all * cover before * add me --------- Co-authored-by: Hendrik Ranocha --- AUTHORS.md | 1 + .../compressible_euler_multicomponent_1d.jl | 160 +++++----- .../compressible_euler_multicomponent_2d.jl | 204 +++++++------ src/equations/compressible_euler_quasi_1d.jl | 52 ++-- test/test_type.jl | 286 +++++++++++++++--- 5 files changed, 458 insertions(+), 245 deletions(-) diff --git a/AUTHORS.md b/AUTHORS.md index 8e3afcf867..54d6321633 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -43,3 +43,4 @@ are listed in alphabetical order: * Michael Schlottke-Lakemper * Toskan Theine * Andrew Winters +* Huiyu Xie diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl index 6338e04c3e..d4f6b421f7 100644 --- a/src/equations/compressible_euler_multicomponent_1d.jl +++ b/src/equations/compressible_euler_multicomponent_1d.jl @@ -123,14 +123,15 @@ A smooth initial condition used for convergence tests in combination with """ function initial_condition_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations1D) + RealT = eltype(x) c = 2 - A = 0.1 + A = convert(RealT, 0.1) L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f ini = c + A * sin(omega * (x[1] - t)) - v1 = 1.0 + v1 = 1 rho = ini @@ -144,7 +145,7 @@ function initial_condition_convergence_test(x, t, prim1 = rho * v1 prim2 = rho^2 - prim_other = SVector{2, real(equations)}(prim1, prim2) + prim_other = SVector(prim1, prim2) return vcat(prim_other, prim_rho) end @@ -159,26 +160,27 @@ Source terms used for convergence tests in combination with @inline function source_terms_convergence_test(u, x, t, equations::CompressibleEulerMulticomponentEquations1D) # Same settings as in `initial_condition` + RealT = eltype(u) c = 2 - A = 0.1 + A = convert(RealT, 0.1) L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f gamma = totalgamma(u, equations) x1, = x si, co = sincos((t - x1) * omega) - tmp = (-((4 * si * A - 4c) + 1) * (gamma - 1) * co * A * omega) / 2 + tmp = (-((4 * si * A - 4 * c) + 1) * (gamma - 1) * co * A * omega) / 2 # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1 - du_rho = SVector{ncomponents(equations), real(equations)}(0.0 + du_rho = SVector{ncomponents(equations), real(equations)}(0 for i in eachcomponent(equations)) du1 = tmp du2 = tmp - du_other = SVector{2, real(equations)}(du1, du2) + du_other = SVector(du1, du2) return vcat(du_other, du_rho) end @@ -194,26 +196,27 @@ A for multicomponent adapted weak blast wave adapted to multicomponent and taken function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerMulticomponentEquations1D) # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) - inicenter = SVector(0.0) + RealT = eltype(x) + inicenter = SVector(0) x_norm = x[1] - inicenter[1] r = abs(x_norm) - cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm) + cos_phi = x_norm > 0 ? 1 : -1 - 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 : + 1 : 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 - 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{2, real(equations)}(v1, p) + prim_other = SVector(v1, p) return prim2cons(vcat(prim_other, prim_rho), equations) end @@ -227,13 +230,13 @@ end v1 = rho_v1 / rho gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * v1^2) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * v1^2) f_rho = densities(u, v1, equations) f1 = rho_v1 * v1 + p f2 = (rho_e + p) * v1 - f_other = SVector{2, real(equations)}(f1, f2) + f_other = SVector(f1, f2) return vcat(f_other, f_rho) end @@ -255,7 +258,7 @@ Entropy conserving two-point flux by rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 2], u_rr[i + 2]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 2] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 2] + u_rr[i + 2]) for i in eachcomponent(equations)) @@ -269,13 +272,14 @@ Entropy conserving two-point flux by # extract velocities v1_ll = rho_v1_ll / rho_ll v1_rr = rho_v1_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) - v1_square = 0.5 * (v1_ll^2 + v1_rr^2) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v1_square = 0.5f0 * (v1_ll^2 + v1_rr^2) v_sum = v1_avg - 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] @@ -283,14 +287,14 @@ Entropy conserving two-point flux by help1_rr += u_rr[i + 2] * cv[i] end - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_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 * (v1_ll^2)) / help1_ll + T_rr = (rho_e_rr - 0.5f0 * rho_rr * (v1_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 - 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)) @@ -299,9 +303,9 @@ Entropy conserving two-point flux by help2 += f_rho[i] end f1 = (help2) * v1_avg + enth / T - f2 = (help1) / T_log - 0.5 * (v1_square) * (help2) + v1_avg * f1 + f2 = (help1) / T_log - 0.5f0 * (v1_square) * (help2) + v1_avg * f1 - f_other = SVector{2, real(equations)}(f1, f2) + f_other = SVector(f1, f2) return vcat(f_other, f_rho) end @@ -330,7 +334,7 @@ See also rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 2], u_rr[i + 2]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 2] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 2] + u_rr[i + 2]) for i in eachcomponent(equations)) @@ -339,25 +343,26 @@ See also rho_rr = density(u_rr, equations) # Calculating gamma - gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + gamma = totalgamma(0.5f0 * (u_ll + u_rr), equations) inv_gamma_minus_one = 1 / (gamma - 1) # extract velocities v1_ll = rho_v1_ll / rho_ll v1_rr = rho_v1_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr) # density flux f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg for i in eachcomponent(equations)) # helpful variables - f_rho_sum = zero(v1_ll) - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) - enth_ll = zero(v1_ll) - enth_rr = zero(v1_rr) + RealT = eltype(u_ll) + f_rho_sum = zero(RealT) + help1_ll = zero(RealT) + help1_rr = zero(RealT) + enth_ll = zero(RealT) + enth_rr = zero(RealT) for i in eachcomponent(equations) enth_ll += u_ll[i + 2] * gas_constants[i] enth_rr += u_rr[i + 2] * gas_constants[i] @@ -367,18 +372,18 @@ See also end # temperature and pressure - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2)) / help1_rr + T_ll = (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2)) / help1_ll + T_rr = (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2)) / help1_rr p_ll = T_ll * enth_ll p_rr = T_rr * enth_rr - p_avg = 0.5 * (p_ll + p_rr) + p_avg = 0.5f0 * (p_ll + p_rr) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) # momentum and energy flux f1 = f_rho_sum * v1_avg + p_avg f2 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (p_ll * v1_rr + p_rr * v1_ll) - f_other = SVector{2, real(equations)}(f1, f2) + 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) + f_other = SVector(f1, f2) return vcat(f_other, f_rho) end @@ -398,8 +403,8 @@ end v_ll = rho_v1_ll / rho_ll v_rr = rho_v1_rr / rho_rr - p_ll = (gamma_ll - 1) * (rho_e_ll - 1 / 2 * rho_ll * v_ll^2) - p_rr = (gamma_rr - 1) * (rho_e_rr - 1 / 2 * rho_rr * v_rr^2) + p_ll = (gamma_ll - 1) * (rho_e_ll - 0.5f0 * rho_ll * v_ll^2) + p_rr = (gamma_rr - 1) * (rho_e_rr - 0.5f0 * rho_rr * v_rr^2) c_ll = sqrt(gamma_ll * p_ll / rho_ll) c_rr = sqrt(gamma_rr * p_rr / rho_rr) @@ -414,7 +419,7 @@ end v1 = rho_v1 / rho gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 1 / 2 * rho * (v1^2)) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2)) c = sqrt(gamma * p / rho) return (abs(v1) + c,) @@ -431,8 +436,8 @@ end v1 = rho_v1 / rho gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2)) - prim_other = SVector{2, real(equations)}(v1, p) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2)) + prim_other = SVector(v1, p) return vcat(prim_other, prim_rho) end @@ -451,7 +456,7 @@ end rho_v1 = rho * v1 - rho_e = p / (gamma - 1) + 0.5 * (rho_v1 * v1) + rho_e = p / (gamma - 1) + 0.5f0 * (rho_v1 * v1) cons_other = SVector{2, RealT}(rho_v1, rho_e) @@ -466,8 +471,9 @@ end rho = density(u, equations) - help1 = zero(rho) - gas_constant = zero(rho) + RealT = eltype(u) + help1 = zero(RealT) + gas_constant = zero(RealT) for i in eachcomponent(equations) help1 += u[i + 2] * cv[i] gas_constant += gas_constants[i] * (u[i + 2] / rho) @@ -477,10 +483,10 @@ end v_square = v1^2 gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * v_square) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * v_square) s = log(p) - gamma * log(rho) - log(gas_constant) rho_p = rho / p - T = (rho_e - 0.5 * rho * v_square) / (help1) + T = (rho_e - 0.5f0 * rho * v_square) / (help1) entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] * (1 - log(T)) + @@ -492,7 +498,7 @@ end w1 = gas_constant * v1 * rho_p w2 = gas_constant * (-rho_p) - entrop_other = SVector{2, real(equations)}(w1, w2) + entrop_other = SVector(w1, w2) return vcat(entrop_other, entrop_rho) end @@ -507,14 +513,15 @@ end (-cv[i] * log(-w[2]) - cp[i] + w[i + 2] - - 0.5 * w[1]^2 / + 0.5f0 * w[1]^2 / w[2])) for i in eachcomponent(equations)) - rho = zero(cons_rho[1]) - help1 = zero(cons_rho[1]) - help2 = zero(cons_rho[1]) - p = zero(cons_rho[1]) + RealT = eltype(w) + rho = zero(RealT) + help1 = zero(RealT) + help2 = zero(RealT) + p = zero(RealT) for i in eachcomponent(equations) rho += cons_rho[i] help1 += cons_rho[i] * cv[i] * gammas[i] @@ -523,8 +530,8 @@ end end u1 = rho * v1 gamma = help1 / help2 - u2 = p / (gamma - 1) + 0.5 * rho * v1^2 - cons_other = SVector{2, real(equations)}(u1, u2) + u2 = p / (gamma - 1) + 0.5f0 * rho * v1^2 + cons_other = SVector(u1, u2) return vcat(cons_other, cons_rho) end @@ -534,7 +541,8 @@ end rho = density(u, equations) T = temperature(u, equations) - total_entropy = zero(u[1]) + RealT = eltype(u) + total_entropy = zero(RealT) for i in eachcomponent(equations) total_entropy -= u[i + 2] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 2])) end @@ -548,7 +556,9 @@ end rho_v1, rho_e = u rho = density(u, equations) - help1 = zero(rho) + + RealT = eltype(u) + help1 = zero(RealT) for i in eachcomponent(equations) help1 += u[i + 2] * cv[i] @@ -556,7 +566,7 @@ end v1 = rho_v1 / rho v_square = v1^2 - T = (rho_e - 0.5 * rho * v_square) / help1 + T = (rho_e - 0.5f0 * rho * v_square) / help1 return T end @@ -570,8 +580,9 @@ partial density fractions as well as the partial specific heats at constant volu @inline function totalgamma(u, equations::CompressibleEulerMulticomponentEquations1D) @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 + 2] * cv[i] * gammas[i] @@ -587,13 +598,14 @@ end rho = density(u, equations) gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * (rho_v1^2) / rho) + p = (gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2) / rho) return p end @inline function density(u, equations::CompressibleEulerMulticomponentEquations1D) - rho = zero(u[1]) + RealT = eltype(u) + rho = zero(RealT) for i in eachcomponent(equations) rho += u[i + 2] diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 60fce222f2..3473f88733 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -127,15 +127,16 @@ A smooth initial condition used for convergence tests in combination with """ function initial_condition_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations2D) + RealT = eltype(x) c = 2 - A = 0.1 + A = convert(RealT, 0.1) L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f ini = c + A * sin(omega * (x[1] + x[2] - t)) - v1 = 1.0 - v2 = 1.0 + v1 = 1 + v2 = 1 rho = ini @@ -150,7 +151,7 @@ function initial_condition_convergence_test(x, t, prim2 = rho * v2 prim3 = rho^2 - prim_other = SVector{3, real(equations)}(prim1, prim2, prim3) + prim_other = SVector(prim1, prim2, prim3) return vcat(prim_other, prim_rho) end @@ -165,11 +166,12 @@ Source terms used for convergence tests in combination with @inline function source_terms_convergence_test(u, x, t, equations::CompressibleEulerMulticomponentEquations2D) # Same settings as in `initial_condition` + RealT = eltype(u) c = 2 - A = 0.1 + A = convert(RealT, 0.1) L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f gamma = totalgamma(u, equations) @@ -191,9 +193,9 @@ Source terms used for convergence tests in combination with du1 = tmp5 du2 = tmp5 - du3 = 2 * ((tmp6 - 1.0) * tmp3 + tmp6 * gamma) * tmp1 + du3 = 2 * ((tmp6 - 1) * tmp3 + tmp6 * gamma) * tmp1 - du_other = SVector{3, real(equations)}(du1, du2, du3) + du_other = SVector(du1, du2, du3) return vcat(du_other, du_rho) end @@ -210,29 +212,30 @@ function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerMulticomponentEquations2D) # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # 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 : + 1 : 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{3, real(equations)}(v1, v2, p) + prim_other = SVector(v1, v2, p) return prim2cons(vcat(prim_other, prim_rho), equations) end @@ -247,7 +250,7 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2 + v2^2)) if orientation == 1 f_rho = densities(u, v1, equations) @@ -261,7 +264,7 @@ end f3 = (rho_e + p) * v2 end - f_other = SVector{3, real(equations)}(f1, f2, f3) + f_other = SVector(f1, f2, f3) return vcat(f_other, f_rho) end @@ -277,14 +280,14 @@ end v2 = rho_v2 / rho v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2 + v2^2)) f_rho = densities(u, v_normal, equations) f1 = rho_v1 * v_normal + p * normal_direction[1] f2 = rho_v2 * v_normal + p * normal_direction[2] f3 = (rho_e + p) * v_normal - f_other = SVector{3, real(equations)}(f1, f2, f3) + f_other = SVector(f1, f2, f3) return vcat(f_other, f_rho) end @@ -306,7 +309,7 @@ Adaption of the entropy conserving two-point flux by rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], u_rr[i + 3]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 3] + u_rr[i + 3]) for i in eachcomponent(equations)) @@ -319,15 +322,16 @@ Adaption of the entropy conserving two-point flux by v2_ll = rho_v2_ll / rho_ll v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v1_square = 0.5 * (v1_ll^2 + v1_rr^2) - v2_square = 0.5 * (v2_ll^2 + v2_rr^2) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v1_square = 0.5f0 * (v1_ll^2 + v1_rr^2) + v2_square = 0.5f0 * (v2_ll^2 + v2_rr^2) v_sum = v1_avg + v2_avg - 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] @@ -335,14 +339,14 @@ Adaption of the entropy conserving two-point flux by help1_rr += u_rr[i + 3] * cv[i] end - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_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 * (v1_ll^2 + v2_ll^2)) / help1_ll + T_rr = (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2 + v2_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 - 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)) @@ -352,7 +356,7 @@ Adaption of the entropy conserving two-point flux by end f1 = (help2) * v1_avg + enth / T f2 = (help2) * v2_avg - f3 = (help1) / T_log - 0.5 * (v1_square + v2_square) * (help2) + v1_avg * f1 + + f3 = (help1) / T_log - 0.5f0 * (v1_square + v2_square) * (help2) + v1_avg * f1 + v2_avg * f2 else f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v2_avg @@ -363,10 +367,10 @@ Adaption of the entropy conserving two-point flux by end f1 = (help2) * v1_avg f2 = (help2) * v2_avg + enth / T - f3 = (help1) / T_log - 0.5 * (v1_square + v2_square) * (help2) + v1_avg * f1 + + f3 = (help1) / T_log - 0.5f0 * (v1_square + v2_square) * (help2) + v1_avg * f1 + v2_avg * f2 end - f_other = SVector{3, real(equations)}(f1, f2, f3) + f_other = SVector(f1, f2, f3) return vcat(f_other, f_rho) end @@ -395,7 +399,7 @@ See also rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], u_rr[i + 3]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 3] + u_rr[i + 3]) for i in eachcomponent(equations)) @@ -404,23 +408,24 @@ See also rho_rr = density(u_rr, equations) # Calculating gamma - gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + gamma = totalgamma(0.5f0 * (u_ll + u_rr), equations) inv_gamma_minus_one = 1 / (gamma - 1) # extract velocities v1_ll = rho_v1_ll / rho_ll v1_rr = rho_v1_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) v2_ll = rho_v2_ll / rho_ll v2_rr = rho_v2_rr / rho_rr - v2_avg = 0.5 * (v2_ll + v2_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) # helpful variables - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) - enth_ll = zero(v1_ll) - enth_rr = zero(v1_rr) + RealT = eltype(u_ll) + help1_ll = zero(RealT) + help1_rr = zero(RealT) + enth_ll = zero(RealT) + enth_rr = zero(RealT) for i in eachcomponent(equations) enth_ll += u_ll[i + 3] * gas_constants[i] enth_rr += u_rr[i + 3] * gas_constants[i] @@ -429,14 +434,14 @@ See also end # temperature and pressure - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr + T_ll = (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll + T_rr = (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr p_ll = T_ll * enth_ll p_rr = T_rr * enth_rr - p_avg = 0.5 * (p_ll + p_rr) + p_avg = 0.5f0 * (p_ll + p_rr) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - f_rho_sum = zero(T_rr) + f_rho_sum = zero(RealT) if orientation == 1 f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg for i in eachcomponent(equations)) @@ -446,7 +451,7 @@ See also f1 = f_rho_sum * v1_avg + p_avg f2 = f_rho_sum * v2_avg f3 = f_rho_sum * (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) else f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v2_avg for i in eachcomponent(equations)) @@ -456,11 +461,11 @@ See also f1 = f_rho_sum * v1_avg f2 = f_rho_sum * v2_avg + p_avg f3 = f_rho_sum * (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) end # momentum and energy flux - f_other = SVector{3, real(equations)}(f1, f2, f3) + f_other = SVector(f1, f2, f3) return vcat(f_other, f_rho) end @@ -474,7 +479,7 @@ end rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], u_rr[i + 3]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 3] + u_rr[i + 3]) for i in eachcomponent(equations)) @@ -483,25 +488,26 @@ end rho_rr = density(u_rr, equations) # Calculating gamma - gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + gamma = totalgamma(0.5f0 * (u_ll + u_rr), equations) inv_gamma_minus_one = 1 / (gamma - 1) # extract velocities v1_ll = rho_v1_ll / rho_ll v1_rr = rho_v1_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) v2_ll = rho_v2_ll / rho_ll v2_rr = rho_v2_rr / rho_rr - v2_avg = 0.5 * (v2_ll + v2_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] # helpful variables - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) - enth_ll = zero(v1_ll) - enth_rr = zero(v1_rr) + RealT = eltype(u_ll) + help1_ll = zero(RealT) + help1_rr = zero(RealT) + enth_ll = zero(RealT) + enth_rr = zero(RealT) for i in eachcomponent(equations) enth_ll += u_ll[i + 3] * gas_constants[i] enth_rr += u_rr[i + 3] * gas_constants[i] @@ -510,15 +516,15 @@ end end # temperature and pressure - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr + T_ll = (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll + T_rr = (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr p_ll = T_ll * enth_ll p_rr = T_rr * enth_rr - p_avg = 0.5 * (p_ll + p_rr) + p_avg = 0.5f0 * (p_ll + p_rr) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - f_rho_sum = zero(T_rr) - f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5 * + f_rho_sum = zero(RealT) + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) for i in eachcomponent(equations)) for i in eachcomponent(equations) @@ -527,7 +533,7 @@ end f1 = f_rho_sum * v1_avg + p_avg * normal_direction[1] f2 = f_rho_sum * v2_avg + p_avg * normal_direction[2] f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) + 0.5f0 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) # momentum and energy flux f_other = SVector(f1, f2, f3) @@ -557,9 +563,9 @@ end end # Compute the sound speeds on the left and right - p_ll = (gamma_ll - 1) * (rho_e_ll - 1 / 2 * (rho_v1_ll^2 + rho_v2_ll^2) / rho_ll) + p_ll = (gamma_ll - 1) * (rho_e_ll - 0.5f0 * (rho_v1_ll^2 + rho_v2_ll^2) / rho_ll) c_ll = sqrt(gamma_ll * p_ll / rho_ll) - p_rr = (gamma_rr - 1) * (rho_e_rr - 1 / 2 * (rho_v1_rr^2 + rho_v2_rr^2) / rho_rr) + p_rr = (gamma_rr - 1) * (rho_e_rr - 0.5f0 * (rho_v1_rr^2 + rho_v2_rr^2) / rho_rr) c_rr = sqrt(gamma_rr * p_rr / rho_rr) λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) @@ -574,7 +580,7 @@ end v2 = rho_v2 / rho gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 1 / 2 * rho * (v1^2 + v2^2)) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2 + v2^2)) c = sqrt(gamma * p / rho) return (abs(v1) + c, abs(v2) + c) @@ -635,8 +641,8 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) - prim_other = SVector{3, real(equations)}(v1, v2, p) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2 + v2^2)) + prim_other = SVector(v1, v2, p) return vcat(prim_other, prim_rho) end @@ -649,8 +655,9 @@ end rho = density(u, equations) # Multicomponent stuff - help1 = zero(rho) - gas_constant = zero(rho) + RealT = eltype(u) + help1 = zero(RealT) + gas_constant = zero(RealT) for i in eachcomponent(equations) help1 += u[i + 3] * cv[i] gas_constant += gas_constants[i] * (u[i + 3] / rho) @@ -661,10 +668,10 @@ end v_square = v1^2 + v2^2 gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * v_square) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * v_square) s = log(p) - gamma * log(rho) - log(gas_constant) rho_p = rho / p - T = (rho_e - 0.5 * rho * v_square) / (help1) + T = (rho_e - 0.5f0 * rho * v_square) / (help1) entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] * (1 - log(T)) + @@ -677,7 +684,7 @@ end w2 = gas_constant * v2 * rho_p w3 = gas_constant * (-rho_p) - entrop_other = SVector{3, real(equations)}(w1, w2, w3) + entrop_other = SVector(w1, w2, w3) return vcat(entrop_other, entrop_rho) end @@ -698,10 +705,11 @@ end 1) for i in eachcomponent(equations)) - rho = zero(cons_rho[1]) - help1 = zero(cons_rho[1]) - help2 = zero(cons_rho[1]) - p = zero(cons_rho[1]) + RealT = eltype(w) + rho = zero(RealT) + help1 = zero(RealT) + help2 = zero(RealT) + p = zero(RealT) for i in eachcomponent(equations) rho += cons_rho[i] help1 += cons_rho[i] * cv[i] * gammas[i] @@ -711,8 +719,8 @@ end u1 = rho * v1 u2 = rho * v2 gamma = help1 / help2 - u3 = p / (gamma - 1) + 0.5 * rho * v_squared - cons_other = SVector{3, real(equations)}(u1, u2, u3) + u3 = p / (gamma - 1) + 0.5f0 * rho * v_squared + cons_other = SVector(u1, u2, u3) return vcat(cons_other, cons_rho) end @@ -728,9 +736,9 @@ end rho_v1 = rho * v1 rho_v2 = rho * v2 - rho_e = p / (gamma - 1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2) + rho_e = p / (gamma - 1) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2) - cons_other = SVector{3, real(equations)}(rho_v1, rho_v2, rho_e) + cons_other = SVector(rho_v1, rho_v2, rho_e) return vcat(cons_other, cons_rho) end @@ -740,7 +748,8 @@ end rho = density(u, equations) T = temperature(u, equations) - total_entropy = zero(u[1]) + RealT = eltype(u) + total_entropy = zero(RealT) for i in eachcomponent(equations) total_entropy -= u[i + 3] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 3])) end @@ -754,7 +763,8 @@ end rho_v1, rho_v2, rho_e = u rho = density(u, equations) - help1 = zero(rho) + RealT = eltype(u) + help1 = zero(RealT) for i in eachcomponent(equations) help1 += u[i + 3] * cv[i] @@ -763,7 +773,7 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v_square = v1^2 + v2^2 - T = (rho_e - 0.5 * rho * v_square) / help1 + T = (rho_e - 0.5f0 * rho * v_square) / help1 return T end @@ -777,8 +787,9 @@ partial density fractions as well as the partial specific heats at constant volu @inline function totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D) @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 + 3] * cv[i] * gammas[i] @@ -794,13 +805,14 @@ end rho = density(u, equations) gamma = totalgamma(u, equations) - rho_times_p = (gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2)) + rho_times_p = (gamma - 1) * (rho * rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2)) return rho_times_p end @inline function density(u, equations::CompressibleEulerMulticomponentEquations2D) - rho = zero(u[1]) + RealT = eltype(u) + rho = zero(RealT) for i in eachcomponent(equations) rho += u[i + 3] diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl index 9c7e3a7269..936487186e 100644 --- a/src/equations/compressible_euler_quasi_1d.jl +++ b/src/equations/compressible_euler_quasi_1d.jl @@ -78,18 +78,19 @@ A smooth initial condition used for convergence tests in combination with """ function initial_condition_convergence_test(x, t, equations::CompressibleEulerEquationsQuasi1D) + RealT = eltype(x) c = 2 - A = 0.1 + A = convert(RealT, 0.1) L = 2 - f = 1 / L - ω = 2 * pi * f + f = 1.0f0 / L + ω = 2 * convert(RealT, pi) * f ini = c + A * sin(ω * (x[1] - t)) rho = ini - v1 = 1.0 + v1 = 1 e = ini^2 / rho - p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2) - a = 1.5 - 0.5 * cos(x[1] * pi) + p = (equations.gamma - 1) * (e - 0.5f0 * rho * v1^2) + a = 1.5f0 - 0.5f0 * cos(x[1] * convert(RealT, pi)) return prim2cons(SVector(rho, v1, p, a), equations) end @@ -108,19 +109,20 @@ as defined in [`initial_condition_convergence_test`](@ref). equations::CompressibleEulerEquationsQuasi1D) # Same settings as in `initial_condition_convergence_test`. # Derivatives calculated with ForwardDiff.jl + RealT = eltype(u) c = 2 - A = 0.1 + A = convert(RealT, 0.1) L = 2 - f = 1 / L - ω = 2 * pi * f + f = 1.0f0 / L + ω = 2 * convert(RealT, pi) * f x1, = x ini(x1, t) = c + A * sin(ω * (x1 - t)) rho(x1, t) = ini(x1, t) - v1(x1, t) = 1.0 + v1(x1, t) = 1 e(x1, t) = ini(x1, t)^2 / rho(x1, t) - p1(x1, t) = (equations.gamma - 1) * (e(x1, t) - 0.5 * rho(x1, t) * v1(x1, t)^2) - a(x1, t) = 1.5 - 0.5 * cos(x1 * pi) + p1(x1, t) = (equations.gamma - 1) * (e(x1, t) - 0.5f0 * rho(x1, t) * v1(x1, t)^2) + a(x1, t) = 1.5f0 - 0.5f0 * cos(x1 * pi) arho(x1, t) = a(x1, t) * rho(x1, t) arhou(x1, t) = arho(x1, t) * v1(x1, t) @@ -142,7 +144,7 @@ as defined in [`initial_condition_convergence_test`](@ref). du2 = darhou_dt(x1, t) + darhouu_dx(x1, t) + a(x1, t) * dp1_dx(x1, t) du3 = daE_dt(x1, t) + dauEp_dx(x1, t) - return SVector(du1, du2, du3, 0.0) + return SVector(du1, du2, du3, 0) end # Calculate 1D flux for a single point @@ -157,7 +159,7 @@ end f2 = a_rho_v1 * v1 f3 = a * v1 * (e + p) - return SVector(f1, f2, f3, zero(eltype(u))) + return SVector(f1, f2, f3, 0) end """ @@ -189,9 +191,7 @@ Further details are available in the paper: # in the arithmetic average of {p}. p_avg = p_ll + p_rr - z = zero(eltype(u_ll)) - - return SVector(z, a_ll * p_avg, z, z) + return SVector(0, a_ll * p_avg, 0, 0) end # While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that @@ -239,19 +239,19 @@ Further details are available in the paper: # 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) - a_v1_avg = 0.5 * (a_ll * v1_ll + a_rr * v1_rr) - p_avg = 0.5 * (p_ll + p_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + a_v1_avg = 0.5f0 * (a_ll * v1_ll + a_rr * v1_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr) # Calculate fluxes # Ignore orientation since it is always "1" in 1D f1 = rho_mean * a_v1_avg f2 = rho_mean * a_v1_avg * v1_avg f3 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll) + 0.5f0 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll) - return SVector(f1, f2, f3, zero(eltype(u_ll))) + return SVector(f1, f2, f3, 0) end # While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that @@ -276,13 +276,13 @@ end e_ll = a_e_ll / a_ll v1_ll = a_rho_v1_ll / a_rho_ll v_mag_ll = abs(v1_ll) - p_ll = (equations.gamma - 1) * (e_ll - 0.5 * rho_ll * v_mag_ll^2) + p_ll = (equations.gamma - 1) * (e_ll - 0.5f0 * rho_ll * v_mag_ll^2) c_ll = sqrt(equations.gamma * p_ll / rho_ll) rho_rr = a_rho_rr / a_rr e_rr = a_e_rr / a_rr v1_rr = a_rho_v1_rr / a_rho_rr v_mag_rr = abs(v1_rr) - p_rr = (equations.gamma - 1) * (e_rr - 0.5 * rho_rr * v_mag_rr^2) + p_rr = (equations.gamma - 1) * (e_rr - 0.5f0 * rho_rr * v_mag_rr^2) c_rr = sqrt(equations.gamma * p_rr / rho_rr) λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr) @@ -293,7 +293,7 @@ end rho = a_rho / a v1 = a_rho_v1 / a_rho e = a_e / a - p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2) + p = (equations.gamma - 1) * (e - 0.5f0 * rho * v1^2) c = sqrt(equations.gamma * p / rho) return (abs(v1) + c,) diff --git a/test/test_type.jl b/test/test_type.jl index de02ec4711..7bea104bb8 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -47,10 +47,13 @@ isdir(outdir) && rm(outdir, recursive = true) for direction in directions if RealT == Float32 # check `surface_flux_function` (test broken) - @test_broken eltype(boundary_condition_wall(u_inner, orientation, - direction, x, t, - surface_flux_function, - equations)) == RealT + @test_broken eltype(@inferred boundary_condition_wall(u_inner, + orientation, + direction, x, + t, + surface_flux_function, + equations)) == + RealT else @test eltype(@inferred boundary_condition_wall(u_inner, orientation, direction, x, t, @@ -62,10 +65,12 @@ isdir(outdir) && rm(outdir, recursive = true) if RealT == Float32 # check `surface_flux_function` (test broken) - @test_broken eltype(boundary_condition_slip_wall(u_inner, normal_direction, - x, t, - surface_flux_function, - equations)) == RealT + @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner, + normal_direction, + x, t, + surface_flux_function, + equations)) == + RealT else @test eltype(@inferred boundary_condition_slip_wall(u_inner, normal_direction, x, t, @@ -75,7 +80,7 @@ isdir(outdir) && rm(outdir, recursive = true) end @test eltype(@inferred flux(u, normal_direction, equations)) == RealT - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == RealT @test eltype(@inferred dissipation(u_ll, u_rr, normal_direction, equations)) == @@ -83,7 +88,7 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations @test eltype(@inferred flux(u, orientation, equations)) == RealT - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @test eltype(@inferred dissipation(u_ll, u_rr, orientation, equations)) == @@ -141,7 +146,8 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred flux_hllc(u_ll, u_rr, orientation, equations)) == RealT if RealT == Float32 # check `ln_mean` (test broken) - @test_broken eltype(flux_chandrashekar(u_ll, u_rr, orientation, equations)) == + @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == RealT else @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, @@ -150,7 +156,8 @@ isdir(outdir) && rm(outdir, recursive = true) end if RealT == Float32 # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(flux_ranocha(u_ll, u_rr, orientation, equations)) == + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, + equations)) == RealT else @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == @@ -167,7 +174,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations))) == RealT - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + @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 @@ -182,11 +189,11 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred prim2cons(u, equations)) == RealT @test eltype(@inferred cons2entropy(u, equations)) == RealT @test eltype(@inferred entropy2cons(u, equations)) == RealT - @test eltype(@inferred density(u, equations)) == RealT - @test eltype(@inferred pressure(u, equations)) == RealT - @test eltype(@inferred density_pressure(u, equations)) == RealT - @test eltype(@inferred entropy(cons, equations)) == RealT - @test eltype(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT end end @@ -247,8 +254,9 @@ isdir(outdir) && rm(outdir, recursive = true) RealT if RealT == Float32 # check `ln_mean` (test broken) - @test_broken eltype(flux_chandrashekar(u_ll, u_rr, normal_direction, - equations)) == + @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, + normal_direction, + equations)) == RealT else @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, normal_direction, @@ -257,7 +265,8 @@ isdir(outdir) && rm(outdir, recursive = true) end if RealT == Float32 # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(flux_ranocha(u_ll, u_rr, normal_direction, equations)) == + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, + equations)) == RealT else @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, @@ -273,7 +282,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations))) == RealT - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == RealT @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, @@ -299,8 +308,9 @@ isdir(outdir) && rm(outdir, recursive = true) RealT if RealT == Float32 # check `ln_mean` (test broken) - @test_broken eltype(flux_chandrashekar(u_ll, u_rr, orientation, - equations)) == + @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, + orientation, + equations)) == RealT else @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, @@ -309,7 +319,8 @@ isdir(outdir) && rm(outdir, recursive = true) end if RealT == Float32 # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(flux_ranocha(u_ll, u_rr, orientation, equations)) == + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, + equations)) == RealT else @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == @@ -329,7 +340,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations))) == RealT - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + @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, @@ -348,15 +359,16 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred prim2cons(u, equations)) == RealT @test eltype(@inferred cons2entropy(u, equations)) == RealT @test eltype(@inferred entropy2cons(u, equations)) == RealT - @test eltype(@inferred Trixi.entropy_guermond_etal(u, equations)) == RealT @test eltype(@inferred Trixi.cons2entropy_guermond_etal(u, equations)) == RealT - @test eltype(@inferred density(u, equations)) == RealT - @test eltype(@inferred pressure(u, equations)) == RealT - @test eltype(@inferred density_pressure(u, equations)) == RealT - @test eltype(@inferred entropy(cons, equations)) == RealT - @test eltype(@inferred Trixi.entropy_math(cons, equations)) == RealT - @test eltype(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT - @test eltype(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_guermond_etal(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_math(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT + # TODO: test `gradient_conservative`, not necessary but good to have end end @@ -416,22 +428,24 @@ isdir(outdir) && rm(outdir, recursive = true) RealT if RealT == Float32 # check `ln_mean` (test broken) - @test_broken eltype(flux_chandrashekar(u_ll, u_rr, normal_direction, - equations)) == RealT + @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, + normal_direction, + equations)) == RealT else @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, normal_direction, equations)) == RealT end if RealT == Float32 # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(flux_ranocha(u_ll, u_rr, normal_direction, equations)) == + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, + equations)) == RealT else @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, equations)) == RealT end - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == RealT @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, @@ -454,15 +468,17 @@ isdir(outdir) && rm(outdir, recursive = true) RealT if RealT == Float32 # check `ln_mean` (test broken) - @test_broken eltype(flux_chandrashekar(u_ll, u_rr, orientation, - equations)) == RealT + @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, + orientation, + equations)) == RealT else @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, equations)) == RealT end if RealT == Float32 # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(flux_ranocha(u_ll, u_rr, orientation, equations)) == + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, + equations)) == RealT else @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == @@ -473,7 +489,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations))) == RealT - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + @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 @@ -488,13 +504,185 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred prim2cons(u, equations)) == RealT @test eltype(@inferred cons2entropy(u, equations)) == RealT @test eltype(@inferred entropy2cons(u, equations)) == RealT - @test eltype(@inferred density(u, equations)) == RealT - @test eltype(@inferred pressure(u, equations)) == RealT - @test eltype(@inferred density_pressure(u, equations)) == RealT - @test eltype(@inferred entropy(cons, equations)) == RealT - @test eltype(@inferred Trixi.entropy_math(cons, equations)) == RealT - @test eltype(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT - @test eltype(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_math(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT + end + end + + @timed_testset "Compressible Euler Multicomponent 1D" begin + for RealT in (Float32, Float64) + gammas = (RealT(1.4), RealT(1.4)) + gas_constants = (RealT(0.4), RealT(0.4)) + equations = @inferred CompressibleEulerMulticomponentEquations1D(gammas = gammas, + gas_constants = gas_constants) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = SVector(one(RealT), one(RealT), one(RealT), one(RealT)) + orientation = 1 + + @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 source_terms_convergence_test(u, 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_chandrashekar(u_ll, u_rr, orientation, + equations)) == RealT + else + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == RealT + end + if RealT == Float32 + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, + equations)) == + RealT + else + @test eltype(@inferred flux_ranocha(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 eltype(@inferred entropy2cons(u, equations)) == RealT + @test typeof(@inferred Trixi.total_entropy(u, equations)) == RealT + @test typeof(@inferred Trixi.temperature(u, equations)) == RealT + @test typeof(@inferred Trixi.totalgamma(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + end + end + + @timed_testset "Compressible Euler Multicomponent 2D" begin + for RealT in (Float32, Float64) + gammas = (RealT(1.4), RealT(1.4)) + gas_constants = (RealT(0.4), RealT(0.4)) + equations = @inferred CompressibleEulerMulticomponentEquations2D(gammas = gammas, + gas_constants = gas_constants) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = SVector(one(RealT), one(RealT), one(RealT), one(RealT), + one(RealT)) + orientations = [1, 2] + normal_direction = SVector(one(RealT), zero(RealT)) + + @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 source_terms_convergence_test(u, x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + if RealT == Float32 + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, + equations)) == RealT + else + @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, + equations)) == RealT + end + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, + orientation, + equations)) == RealT + else + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == RealT + end + if RealT == Float32 + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, + equations)) == RealT + else + @test eltype(@inferred flux_ranocha(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 eltype(@inferred entropy2cons(u, equations)) == RealT + @test typeof(@inferred Trixi.total_entropy(u, equations)) == RealT + @test typeof(@inferred Trixi.temperature(u, equations)) == RealT + @test typeof(@inferred Trixi.totalgamma(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + end + end + + @timed_testset "Compressible Euler Quasi 1D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleEulerEquationsQuasi1D(RealT(1.4)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = SVector(one(RealT), one(RealT), one(RealT), one(RealT)) + orientation = 1 + normal_direction = normal_ll = normal_rr = SVector(one(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_nonconservative_chan_etal(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_nonconservative_chan_etal(u_ll, u_rr, + normal_direction, + equations)) == + RealT + @test eltype(@inferred flux_nonconservative_chan_etal(u_ll, u_rr, normal_ll, + normal_rr, equations)) == + RealT + if RealT == Float32 + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_chan_etal(u_ll, u_rr, orientation, + equations)) == + RealT + else + @test eltype(@inferred flux_chan_etal(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 entropy(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT end end end