From 866530092e95332b2fd51a7063d3140ba72a577c Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Mon, 3 Jun 2024 05:19:32 -0700 Subject: [PATCH] Add numerical support of other real types (`Float32`) (#1909) * start with src/equations * revise after the first review * format src/equations * another version * revise after the second review * apply suggestions from code review - comments Co-authored-by: Hendrik Ranocha * remove TODO * fix small errors * complete compressible_euler * fix error and format * revise min max * revise based on new docs * revise based on new comments * start sample test * add more tests * fix typos and comments * change code review * add to CI * apply inferred macro Co-authored-by: Hendrik Ranocha * complete * add Trixi --------- Co-authored-by: Hendrik Ranocha --- src/equations/acoustic_perturbation_2d.jl | 56 +-- src/equations/compressible_euler_1d.jl | 213 ++++----- src/equations/compressible_euler_2d.jl | 450 +++++++++---------- src/equations/compressible_euler_3d.jl | 379 ++++++++-------- test/runtests.jl | 1 + test/test_type.jl | 502 ++++++++++++++++++++++ 6 files changed, 1075 insertions(+), 526 deletions(-) create mode 100644 test/test_type.jl diff --git a/src/equations/acoustic_perturbation_2d.jl b/src/equations/acoustic_perturbation_2d.jl index f4ce770e1e9..1bde94a648a 100644 --- a/src/equations/acoustic_perturbation_2d.jl +++ b/src/equations/acoustic_perturbation_2d.jl @@ -112,9 +112,9 @@ A constant initial condition where the state variables are zero and the mean flo Uses the global mean values from `equations`. """ function initial_condition_constant(x, t, equations::AcousticPerturbationEquations2D) - v1_prime = 0.0 - v2_prime = 0.0 - p_prime_scaled = 0.0 + v1_prime = 0 + v2_prime = 0 + p_prime_scaled = 0 return SVector(v1_prime, v2_prime, p_prime_scaled, global_mean_vars(equations)...) end @@ -127,12 +127,13 @@ A smooth initial condition used for convergence tests in combination with """ function initial_condition_convergence_test(x, t, equations::AcousticPerturbationEquations2D) - c = 2.0 - A = 0.2 - L = 2.0 - f = 2.0 / L - a = 1.0 - omega = 2 * pi * f + RealT = eltype(x) + a = 1 + c = 2 + L = 2 + f = 2.0f0 / L + A = convert(RealT, 0.2) + omega = 2 * convert(RealT, pi) * f init = c + A * sin(omega * (x[1] + x[2] - a * t)) v1_prime = init @@ -154,12 +155,13 @@ function source_terms_convergence_test(u, x, t, equations::AcousticPerturbationEquations2D) v1_mean, v2_mean, c_mean, rho_mean = cons2mean(u, equations) - c = 2.0 - A = 0.2 - L = 2.0 - f = 2.0 / L - a = 1.0 - omega = 2 * pi * f + RealT = eltype(u) + a = 1 + c = 2 + L = 2 + f = 2.0f0 / L + A = convert(RealT, 0.2) + omega = 2 * convert(RealT, pi) * f si, co = sincos(omega * (x[1] + x[2] - a * t)) tmp = v1_mean + v2_mean - a @@ -168,7 +170,7 @@ function source_terms_convergence_test(u, x, t, du3 = A * omega * co * (2 * c_mean^2 * rho_mean + 2 * c * tmp + 2 * A * tmp * si) / c_mean^2 - du4 = du5 = du6 = du7 = 0.0 + du4 = du5 = du6 = du7 = 0 return SVector(du1, du2, du3, du4, du5, du6, du7) end @@ -179,8 +181,8 @@ end A Gaussian pulse in a constant mean flow. Uses the global mean values from `equations`. """ function initial_condition_gauss(x, t, equations::AcousticPerturbationEquations2D) - v1_prime = 0.0 - v2_prime = 0.0 + v1_prime = 0 + v2_prime = 0 p_prime = exp(-4 * (x[1]^2 + x[2]^2)) prim = SVector(v1_prime, v2_prime, p_prime, global_mean_vars(equations)...) @@ -240,8 +242,8 @@ function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, u_normal = normal[1] * u_inner[1] + normal[2] * u_inner[2] # create the "external" boundary solution state - u_boundary = SVector(u_inner[1] - 2.0 * u_normal * normal[1], - u_inner[2] - 2.0 * u_normal * normal[2], + u_boundary = SVector(u_inner[1] - 2 * u_normal * normal[1], + u_inner[2] - 2 * u_normal * normal[2], u_inner[3], cons2mean(u_inner, equations)...) # calculate the boundary flux @@ -257,13 +259,14 @@ end v1_mean, v2_mean, c_mean, rho_mean = cons2mean(u, equations) # Calculate flux for conservative state variables + RealT = eltype(u) if orientation == 1 f1 = v1_mean * v1_prime + v2_mean * v2_prime + c_mean^2 * p_prime_scaled / rho_mean - f2 = zero(eltype(u)) + f2 = zero(RealT) f3 = rho_mean * v1_prime + v1_mean * p_prime_scaled else - f1 = zero(eltype(u)) + f1 = zero(RealT) f2 = v1_mean * v1_prime + v2_mean * v2_prime + c_mean^2 * p_prime_scaled / rho_mean f3 = rho_mean * v2_prime + v2_mean * p_prime_scaled @@ -272,7 +275,7 @@ end # The rest of the state variables are actually variable coefficients, hence the flux should be # zero. See https://github.com/trixi-framework/Trixi.jl/issues/358#issuecomment-784828762 # for details. - f4 = f5 = f6 = f7 = zero(eltype(u)) + f4 = f5 = f6 = f7 = 0 return SVector(f1, f2, f3, f4, f5, f6, f7) end @@ -313,7 +316,7 @@ end # The rest of the state variables are actually variable coefficients, hence the flux should be # zero. See https://github.com/trixi-framework/Trixi.jl/issues/358#issuecomment-784828762 # for details. - f4 = f5 = f6 = f7 = zero(eltype(u)) + f4 = f5 = f6 = f7 = 0 return SVector(f1, f2, f3, f4, f5, f6, f7) end @@ -344,8 +347,9 @@ end equations::AcousticPerturbationEquations2D) λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations) - diss = -0.5 * λ * (u_rr - u_ll) - z = zero(eltype(u_ll)) + diss = -0.5f0 * λ * (u_rr - u_ll) + z = 0 + return SVector(diss[1], diss[2], diss[3], z, z, z, z) end diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index a50c896cd1c..a71750ff98c 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -53,9 +53,10 @@ varnames(::typeof(cons2prim), ::CompressibleEulerEquations1D) = ("rho", "v1", "p A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equations::CompressibleEulerEquations1D) - rho = 1.0 - rho_v1 = 0.1 - rho_e = 10.0 + RealT = eltype(x) + rho = 1 + rho_v1 = convert(RealT, 0.1) + rho_e = 10 return SVector(rho, rho_v1, rho_e) end @@ -68,11 +69,12 @@ A smooth initial condition used for convergence tests in combination with """ function initial_condition_convergence_test(x, t, equations::CompressibleEulerEquations1D) + 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 @@ -92,11 +94,12 @@ Source terms used for convergence tests in combination with @inline function source_terms_convergence_test(u, x, t, equations::CompressibleEulerEquations1D) # Same settings as in `initial_condition` + 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 γ = equations.gamma x1, = x @@ -108,8 +111,8 @@ Source terms used for convergence tests in combination with # Note that d/dt rho = -d/dx rho. # This yields du2 = du3 = d/dx p (derivative of pressure). # Other terms vanish because of v = 1. - du1 = zero(eltype(u)) - du2 = rho_x * (2 * rho - 0.5) * (γ - 1) + du1 = 0 + du2 = rho_x * (2 * rho - 0.5f0) * (γ - 1) du3 = du2 return SVector(du1, du2, du3) @@ -130,11 +133,12 @@ with the following parameters - polydeg = 5 """ function initial_condition_density_wave(x, t, equations::CompressibleEulerEquations1D) - v1 = 0.1 - rho = 1 + 0.98 * sinpi(2 * (x[1] - t * v1)) + RealT = eltype(x) + v1 = convert(RealT, 0.1) + rho = 1 + convert(RealT, 0.98) * sinpi(2 * (x[1] - t * v1)) rho_v1 = rho * v1 p = 20 - rho_e = p / (equations.gamma - 1) + 1 / 2 * rho * v1^2 + rho_e = p / (equations.gamma - 1) + 0.5f0 * rho * v1^2 return SVector(rho, rho_v1, rho_e) end @@ -150,19 +154,20 @@ function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations1D) # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # Set up polar coordinates - inicenter = SVector(0.0) + RealT = eltype(x) + inicenter = SVector(0) x_norm = x[1] - inicenter[1] r = abs(x_norm) # The following code is equivalent to # phi = atan(0.0, x_norm) # cos_phi = cos(phi) # in 1D but faster - cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm) + cos_phi = x_norm > 0 ? 1 : -1 # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - p = r > 0.5 ? 1.0 : 1.245 + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) return prim2cons(SVector(rho, v1, p), equations) end @@ -183,17 +188,18 @@ with self-gravity from function initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::CompressibleEulerEquations1D) # OBS! this assumes that γ = 2 other manufactured source terms are incorrect - if equations.gamma != 2.0 + if equations.gamma != 2 error("adiabatic constant must be 2 for the coupling convergence test") end - c = 2.0 - A = 0.1 + RealT = eltype(x) + c = 2 + A = convert(RealT, 0.1) ini = c + A * sinpi(x[1] - t) - G = 1.0 # gravitational constant + G = 1 # gravitational constant rho = ini - v1 = 1.0 - p = 2 * ini^2 * G / pi # * 2 / ndims, but ndims==1 here + v1 = 1 + p = 2 * ini^2 * G / convert(RealT, pi) # * 2 / ndims, but ndims==1 here return prim2cons(SVector(rho, v1, p), equations) end @@ -229,31 +235,29 @@ are available in the paper: # Eleuterio F. Toro (2009) # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) - if v_normal <= 0.0 + if v_normal <= 0 sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed p_star = p_local * - (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * - equations.gamma * - equations.inv_gamma_minus_one) - else # v_normal > 0.0 + (1 + 0.5f0 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * + equations.gamma * + equations.inv_gamma_minus_one) + else # v_normal > 0 A = 2 / ((equations.gamma + 1) * rho_local) B = p_local * (equations.gamma - 1) / (equations.gamma + 1) p_star = p_local + - 0.5 * v_normal / A * + 0.5f0 * v_normal / A * (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) end # For the slip wall we directly set the flux as the normal velocity is zero - return SVector(zero(eltype(u_inner)), - p_star, - zero(eltype(u_inner))) + return SVector(0, p_star, 0) end # Calculate 1D flux for a single point @inline function flux(u, orientation::Integer, equations::CompressibleEulerEquations1D) rho, rho_v1, rho_e = u v1 = rho_v1 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho_v1 * v1) # Ignore orientation since it is always "1" in 1D f1 = rho_v1 f2 = rho_v1 * v1 + p @@ -283,14 +287,14 @@ The modification is in the energy flux to guarantee pressure equilibrium and was rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations) # Average each factor of products in flux - rho_avg = 1 / 2 * (rho_ll + rho_rr) - v1_avg = 1 / 2 * (v1_ll + v1_rr) - p_avg = 1 / 2 * (p_ll + p_rr) - kin_avg = 1 / 2 * (v1_ll * v1_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + kin_avg = 0.5f0 * (v1_ll * v1_rr) # Calculate fluxes # Ignore orientation since it is always "1" in 1D - pv1_avg = 1 / 2 * (p_ll * v1_rr + p_rr * v1_ll) + pv1_avg = 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) f1 = rho_avg * v1_avg f2 = f1 * v1_avg + p_avg f3 = p_avg * v1_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv1_avg @@ -316,10 +320,10 @@ Kinetic energy preserving two-point flux by rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations) # Average each factor of products in flux - rho_avg = 1 / 2 * (rho_ll + rho_rr) - v1_avg = 1 / 2 * (v1_ll + v1_rr) - p_avg = 1 / 2 * (p_ll + p_rr) - e_avg = 1 / 2 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + e_avg = 0.5f0 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) # Ignore orientation since it is always "1" in 1D f1 = rho_avg * v1_avg @@ -343,25 +347,25 @@ Entropy conserving two-point flux by # Unpack left and right state rho_ll, v1_ll, p_ll = cons2prim(u_ll, equations) rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations) - beta_ll = 0.5 * rho_ll / p_ll - beta_rr = 0.5 * rho_rr / p_rr - specific_kin_ll = 0.5 * (v1_ll^2) - specific_kin_rr = 0.5 * (v1_rr^2) + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr + specific_kin_ll = 0.5f0 * (v1_ll^2) + specific_kin_rr = 0.5f0 * (v1_rr^2) # Compute the necessary mean values - rho_avg = 0.5 * (rho_ll + rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) rho_mean = ln_mean(rho_ll, rho_rr) beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - p_mean = 0.5 * rho_avg / beta_avg + beta_avg = 0.5f0 * (beta_ll + beta_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + p_mean = 0.5f0 * rho_avg / beta_avg velocity_square_avg = specific_kin_ll + specific_kin_rr # Calculate fluxes # Ignore orientation since it is always "1" in 1D f1 = rho_mean * v1_avg f2 = f1 * v1_avg + p_mean - f3 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f3 = f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + f2 * v1_avg return SVector(f1, f2, f3) @@ -394,16 +398,16 @@ See also # 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) - 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) + 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 * v1_avg f2 = f1 * v1_avg + p_avg f3 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (p_ll * v1_rr + p_rr * v1_ll) + 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) return SVector(f1, f2, f3) end @@ -453,7 +457,7 @@ end equations::CompressibleEulerEquations1D) rho, rho_v1, rho_e = u v1 = rho_v1 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho_v1 * v1) a = sqrt(equations.gamma * p / rho) lambda1 = v1 @@ -466,10 +470,10 @@ end alpha_p = 2 * (equations.gamma - 1) * lambda1_p + lambda2_p + lambda3_p - rho_2gamma = 0.5 * rho / equations.gamma + rho_2gamma = 0.5f0 * rho / equations.gamma f1p = rho_2gamma * alpha_p f2p = rho_2gamma * (alpha_p * v1 + a * (lambda2_p - lambda3_p)) - f3p = rho_2gamma * (alpha_p * 0.5 * v1^2 + a * v1 * (lambda2_p - lambda3_p) + f3p = rho_2gamma * (alpha_p * 0.5f0 * v1^2 + a * v1 * (lambda2_p - lambda3_p) + a^2 * (lambda2_p + lambda3_p) * equations.inv_gamma_minus_one) return SVector(f1p, f2p, f3p) @@ -479,7 +483,7 @@ end equations::CompressibleEulerEquations1D) rho, rho_v1, rho_e = u v1 = rho_v1 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho_v1 * v1) a = sqrt(equations.gamma * p / rho) lambda1 = v1 @@ -492,10 +496,10 @@ end alpha_m = 2 * (equations.gamma - 1) * lambda1_m + lambda2_m + lambda3_m - rho_2gamma = 0.5 * rho / equations.gamma + rho_2gamma = 0.5f0 * rho / equations.gamma f1m = rho_2gamma * alpha_m f2m = rho_2gamma * (alpha_m * v1 + a * (lambda2_m - lambda3_m)) - f3m = rho_2gamma * (alpha_m * 0.5 * v1^2 + a * v1 * (lambda2_m - lambda3_m) + f3m = rho_2gamma * (alpha_m * 0.5f0 * v1^2 + a * v1 * (lambda2_m - lambda3_m) + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one) return SVector(f1m, f2m, f3m) @@ -546,7 +550,7 @@ end equations::CompressibleEulerEquations1D) rho, rho_v1, rho_e = u v1 = rho_v1 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho_v1 * v1) # sound speed and enthalpy a = sqrt(equations.gamma * p / rho) @@ -555,9 +559,9 @@ end # signed Mach number M = v1 / a - p_plus = 0.5 * (1 + equations.gamma * M) * p + p_plus = 0.5f0 * (1 + equations.gamma * M) * p - f1p = 0.25 * rho * a * (M + 1)^2 + f1p = 0.25f0 * rho * a * (M + 1)^2 f2p = f1p * v1 + p_plus f3p = f1p * H @@ -568,7 +572,7 @@ end equations::CompressibleEulerEquations1D) rho, rho_v1, rho_e = u v1 = rho_v1 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho_v1 * v1) # sound speed and enthalpy a = sqrt(equations.gamma * p / rho) @@ -577,9 +581,9 @@ end # signed Mach number M = v1 / a - p_minus = 0.5 * (1 - equations.gamma * M) * p + p_minus = 0.5f0 * (1 - equations.gamma * M) * p - f1m = -0.25 * rho * a * (M - 1)^2 + f1m = -0.25f0 * rho * a * (M - 1)^2 f2m = f1m * v1 + p_minus f3m = f1m * H @@ -634,7 +638,7 @@ end equations::CompressibleEulerEquations1D) rho, rho_v1, rho_e = u v1 = rho_v1 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho_v1 * v1) # sound speed and enthalpy a = sqrt(equations.gamma * p / rho) @@ -644,13 +648,13 @@ end M = v1 / a P = 2 - mu = 1.0 - nu = 0.75 - omega = 2.0 # adjusted from suggested value of 1.5 + mu = 1 + nu = 0.75f0 + omega = 2 # adjusted from suggested value of 1.5 - p_plus = 0.25 * ((M + 1)^2 * (2 - M) - nu * M * (M^2 - 1)^P) * p + p_plus = 0.25f0 * ((M + 1)^2 * (2 - M) - nu * M * (M^2 - 1)^P) * p - f1p = 0.25 * rho * a * ((M + 1)^2 - mu * (M^2 - 1)^P) + f1p = 0.25f0 * rho * a * ((M + 1)^2 - mu * (M^2 - 1)^P) f2p = f1p * v1 + p_plus f3p = f1p * H - omega * rho * a^3 * M^2 * (M^2 - 1)^2 @@ -661,7 +665,7 @@ end equations::CompressibleEulerEquations1D) rho, rho_v1, rho_e = u v1 = rho_v1 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho_v1 * v1) # sound speed and enthalpy a = sqrt(equations.gamma * p / rho) @@ -671,13 +675,13 @@ end M = v1 / a P = 2 - mu = 1.0 - nu = 0.75 - omega = 2.0 # adjusted from suggested value of 1.5 + mu = 1 + nu = 0.75f0 + omega = 2 # adjusted from suggested value of 1.5 - p_minus = 0.25 * ((M - 1)^2 * (2 + M) + nu * M * (M^2 - 1)^P) * p + p_minus = 0.25f0 * ((M - 1)^2 * (2 + M) + nu * M * (M^2 - 1)^P) * p - f1m = -0.25 * rho * a * ((M - 1)^2 - mu * (M^2 - 1)^P) + f1m = -0.25f0 * rho * a * ((M - 1)^2 - mu * (M^2 - 1)^P) f2m = f1m * v1 + p_minus f3m = f1m * H + omega * rho * a^3 * M^2 * (M^2 - 1)^2 @@ -694,11 +698,11 @@ end # Calculate primitive variables and speed of sound v1_ll = rho_v1_ll / rho_ll v_mag_ll = abs(v1_ll) - p_ll = (equations.gamma - 1) * (rho_e_ll - 1 / 2 * rho_ll * v_mag_ll^2) + p_ll = (equations.gamma - 1) * (rho_e_ll - 0.5f0 * rho_ll * v_mag_ll^2) c_ll = sqrt(equations.gamma * p_ll / rho_ll) v1_rr = rho_v1_rr / rho_rr v_mag_rr = abs(v1_rr) - p_rr = (equations.gamma - 1) * (rho_e_rr - 1 / 2 * rho_rr * v_mag_rr^2) + p_rr = (equations.gamma - 1) * (rho_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) @@ -746,12 +750,12 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, v1_ll = rho_v1_ll / rho_ll e_ll = rho_e_ll / rho_ll - p_ll = (equations.gamma - 1) * (rho_e_ll - 1 / 2 * rho_ll * v1_ll^2) + p_ll = (equations.gamma - 1) * (rho_e_ll - 0.5f0 * rho_ll * v1_ll^2) c_ll = sqrt(equations.gamma * p_ll / rho_ll) v1_rr = rho_v1_rr / rho_rr e_rr = rho_e_rr / rho_rr - p_rr = (equations.gamma - 1) * (rho_e_rr - 1 / 2 * rho_rr * v1_rr^2) + p_rr = (equations.gamma - 1) * (rho_e_rr - 0.5f0 * rho_rr * v1_rr^2) c_rr = sqrt(equations.gamma * p_rr / rho_rr) # Obtain left and right fluxes @@ -765,7 +769,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, vel_L = v1_ll vel_R = v1_rr vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho - ekin_roe = 0.5 * vel_roe^2 + ekin_roe = 0.5f0 * vel_roe^2 H_ll = (rho_e_ll + p_ll) / rho_ll H_rr = (rho_e_rr + p_rr) / rho_rr H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho @@ -775,18 +779,18 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, Ssr = max(vel_R + c_rr, vel_roe + c_roe) sMu_L = Ssl - vel_L sMu_R = Ssr - vel_R - if Ssl >= 0.0 + if Ssl >= 0 f1 = f_ll[1] f2 = f_ll[2] f3 = f_ll[3] - elseif Ssr <= 0.0 + elseif Ssr <= 0 f1 = f_rr[1] f2 = f_rr[2] f3 = f_rr[3] else SStar = (p_rr - p_ll + rho_ll * vel_L * sMu_L - rho_rr * vel_R * sMu_R) / (rho_ll * sMu_L - rho_rr * sMu_R) - if Ssl <= 0.0 <= SStar + if Ssl <= 0 <= SStar densStar = rho_ll * sMu_L / (Ssl - SStar) enerStar = e_ll + (SStar - vel_L) * (SStar + p_ll / (rho_ll * sMu_L)) UStar1 = densStar @@ -853,15 +857,15 @@ Compactly summarized: v_roe_mag = v_roe^2 H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) * inv_sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag)) + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5f0 * v_roe_mag)) # Compute convenience constant for positivity preservation, see # https://doi.org/10.1016/0021-9991(91)90211-3 - beta = sqrt(0.5 * (equations.gamma - 1) / equations.gamma) + beta = sqrt(0.5f0 * (equations.gamma - 1) / equations.gamma) # Estimate the edges of the Riemann fan (with positivity conservation) - SsL = min(v_roe - c_roe, v_ll - beta * c_ll, zero(v_roe)) - SsR = max(v_roe + c_roe, v_rr + beta * c_rr, zero(v_roe)) + SsL = min(v_roe - c_roe, v_ll - beta * c_ll, 0) + SsR = max(v_roe + c_roe, v_rr + beta * c_rr, 0) return SsL, SsR end @@ -869,7 +873,7 @@ end @inline function max_abs_speeds(u, equations::CompressibleEulerEquations1D) rho, rho_v1, rho_e = u v1 = rho_v1 / rho - p = (equations.gamma - 1) * (rho_e - 1 / 2 * rho * v1^2) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho * v1^2) c = sqrt(equations.gamma * p / rho) return (abs(v1) + c,) @@ -880,7 +884,7 @@ end rho, rho_v1, rho_e = u v1 = rho_v1 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1 * v1) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho_v1 * v1) return SVector(rho, v1, p) end @@ -891,11 +895,12 @@ end v1 = rho_v1 / rho v_square = v1^2 - p = (equations.gamma - 1) * (rho_e - 0.5 * rho * v_square) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho * v_square) s = log(p) - equations.gamma * log(rho) rho_p = rho / p - w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - 0.5 * rho_p * v_square + w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - + 0.5f0 * rho_p * v_square w2 = rho_p * v1 w3 = -rho_p @@ -912,7 +917,7 @@ end V1, V2, V5 = w .* (gamma - 1) # specific entropy, eq. (53) - s = gamma - V1 + 0.5 * (V2^2) / V5 + s = gamma - V1 + 0.5f0 * (V2^2) / V5 # eq. (52) energy_internal = ((gamma - 1) / (-V5)^gamma)^(equations.inv_gamma_minus_one) * @@ -921,7 +926,7 @@ end # eq. (51) rho = -V5 * energy_internal rho_v1 = V2 * energy_internal - rho_e = (1 - 0.5 * (V2^2) / V5) * energy_internal + rho_e = (1 - 0.5f0 * (V2^2) / V5) * energy_internal return SVector(rho, rho_v1, rho_e) end @@ -929,7 +934,7 @@ end @inline function prim2cons(prim, equations::CompressibleEulerEquations1D) rho, v1, p = prim rho_v1 = rho * v1 - rho_e = p * equations.inv_gamma_minus_one + 0.5 * (rho_v1 * v1) + rho_e = p * equations.inv_gamma_minus_one + 0.5f0 * (rho_v1 * v1) return SVector(rho, rho_v1, rho_e) end @@ -940,20 +945,20 @@ end @inline function pressure(u, equations::CompressibleEulerEquations1D) rho, rho_v1, rho_e = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2) / rho) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2) / rho) return p end @inline function density_pressure(u, equations::CompressibleEulerEquations1D) rho, rho_v1, rho_e = u - rho_times_p = (equations.gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2)) + rho_times_p = (equations.gamma - 1) * (rho * rho_e - 0.5f0 * (rho_v1^2)) return rho_times_p end # Calculate thermodynamic entropy for a conservative state `cons` @inline function entropy_thermodynamic(cons, equations::CompressibleEulerEquations1D) # Pressure - p = (equations.gamma - 1) * (cons[3] - 1 / 2 * (cons[2]^2) / cons[1]) + p = (equations.gamma - 1) * (cons[3] - 0.5f0 * (cons[2]^2) / cons[1]) # Thermodynamic entropy s = log(p) - equations.gamma * log(cons[1]) @@ -980,7 +985,7 @@ end # Calculate kinetic energy for a conservative state `cons` @inline function energy_kinetic(cons, equations::CompressibleEulerEquations1D) - return 0.5 * (cons[2]^2) / cons[1] + return 0.5f0 * (cons[2]^2) / cons[1] end # Calculate internal energy for a conservative state `cons` diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index d15c5c65355..e964d9c2b21 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -60,10 +60,11 @@ varnames(::typeof(cons2prim), ::CompressibleEulerEquations2D) = ("rho", "v1", "v A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equations::CompressibleEulerEquations2D) - rho = 1.0 - rho_v1 = 0.1 - rho_v2 = -0.2 - rho_e = 10.0 + RealT = eltype(x) + rho = 1 + rho_v1 = convert(RealT, 0.1) + rho_v2 = convert(RealT, -0.2) + rho_e = 10 return SVector(rho, rho_v1, rho_v2, rho_e) end @@ -76,11 +77,12 @@ A smooth initial condition used for convergence tests in combination with """ function initial_condition_convergence_test(x, t, equations::CompressibleEulerEquations2D) + 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] + x[2] - t)) rho = ini @@ -101,11 +103,12 @@ Source terms used for convergence tests in combination with @inline function source_terms_convergence_test(u, x, t, equations::CompressibleEulerEquations2D) # Same settings as in `initial_condition` + 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 γ = equations.gamma x1, x2 = x @@ -139,13 +142,14 @@ with the following parameters - polydeg = 5 """ function initial_condition_density_wave(x, t, equations::CompressibleEulerEquations2D) - v1 = 0.1 - v2 = 0.2 - rho = 1 + 0.98 * sinpi(2 * (x[1] + x[2] - t * (v1 + v2))) + RealT = eltype(x) + v1 = convert(RealT, 0.1) + v2 = convert(RealT, 0.2) + rho = 1 + convert(RealT, 0.98) * sinpi(2 * (x[1] + x[2] - t * (v1 + v2))) rho_v1 = rho * v1 rho_v2 = rho * v2 p = 20 - rho_e = p / (equations.gamma - 1) + 1 / 2 * rho * (v1^2 + v2^2) + rho_e = p / (equations.gamma - 1) + 0.5f0 * rho * (v1^2 + v2^2) return SVector(rho, rho_v1, rho_v2, rho_e) end @@ -161,7 +165,7 @@ function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations2D) # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # Set up polar coordinates - inicenter = SVector(0.0, 0.0) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) @@ -169,10 +173,11 @@ function initial_condition_weak_blast_wave(x, t, sin_phi, cos_phi = sincos(phi) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - 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 + RealT = eltype(x) + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + 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) return prim2cons(SVector(rho, v1, v2, p), equations) end @@ -190,18 +195,19 @@ or [`source_terms_eoc_test_euler`](@ref). function initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::CompressibleEulerEquations2D) # OBS! this assumes that γ = 2 other manufactured source terms are incorrect - if equations.gamma != 2.0 + if equations.gamma != 2 error("adiabatic constant must be 2 for the coupling convergence test") end - c = 2.0 - A = 0.1 - ini = c + A * sin(pi * (x[1] + x[2] - t)) - G = 1.0 # gravitational constant + RealT = eltype(x) + c = 2 + A = convert(RealT, 0.1) + ini = c + A * sin(convert(RealT, pi) * (x[1] + x[2] - t)) + G = 1 # gravitational constant rho = ini - v1 = 1.0 - v2 = 1.0 - p = ini^2 * G / pi # * 2 / ndims, but ndims==2 here + v1 = 1 + v2 = 1 + p = ini^2 * G / convert(RealT, pi) # * 2 / ndims, but ndims==2 here return prim2cons(SVector(rho, v1, v2, p), equations) end @@ -218,20 +224,21 @@ in combination with [`initial_condition_eoc_test_coupled_euler_gravity`](@ref). @inline function source_terms_eoc_test_coupled_euler_gravity(u, x, t, equations::CompressibleEulerEquations2D) # Same settings as in `initial_condition_eoc_test_coupled_euler_gravity` - c = 2.0 - A = 0.1 - G = 1.0 # gravitational constant, must match coupling solver - C_grav = -2 * G / pi # 2 == 4 / ndims + RealT = eltype(u) + c = 2 + A = convert(RealT, 0.1) + G = 1 # gravitational constant, must match coupling solver + C_grav = -2 * G / convert(RealT, pi) # 2 == 4 / ndims x1, x2 = x - si, co = sincos(pi * (x1 + x2 - t)) - rhox = A * pi * co + si, co = sincos(convert(RealT, pi) * (x1 + x2 - t)) + rhox = A * convert(RealT, pi) * co rho = c + A * si du1 = rhox du2 = rhox du3 = rhox - du4 = (1.0 - C_grav * rho) * rhox + du4 = (1 - C_grav * rho) * rhox return SVector(du1, du2, du3, du4) end @@ -248,14 +255,15 @@ in combination with [`initial_condition_eoc_test_coupled_euler_gravity`](@ref). @inline function source_terms_eoc_test_euler(u, x, t, equations::CompressibleEulerEquations2D) # Same settings as in `initial_condition_eoc_test_coupled_euler_gravity` - c = 2.0 - A = 0.1 - G = 1.0 - C_grav = -2 * G / pi # 2 == 4 / ndims + RealT = eltype(u) + c = 2 + A = convert(RealT, 0.1) + G = 1 + C_grav = -2 * G / convert(RealT, pi) # 2 == 4 / ndims x1, x2 = x - si, co = sincos(pi * (x1 + x2 - t)) - rhox = A * pi * co + si, co = sincos(convert(RealT, pi) * (x1 + x2 - t)) + rhox = A * convert(RealT, pi) * co rho = c + A * si du1 = rhox @@ -307,25 +315,25 @@ Should be used together with [`UnstructuredMesh2D`](@ref). # Eleuterio F. Toro (2009) # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) - if v_normal <= 0.0 + if v_normal <= 0 sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed p_star = p_local * - (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * - equations.gamma * - equations.inv_gamma_minus_one) - else # v_normal > 0.0 + (1 + 0.5f0 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * + equations.gamma * + equations.inv_gamma_minus_one) + else # v_normal > 0 A = 2 / ((equations.gamma + 1) * rho_local) B = p_local * (equations.gamma - 1) / (equations.gamma + 1) p_star = p_local + - 0.5 * v_normal / A * + 0.5f0 * v_normal / A * (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) end # For the slip wall we directly set the flux as the normal velocity is zero - return SVector(zero(eltype(u_inner)), + return SVector(0, p_star * normal[1], p_star * normal[2], - zero(eltype(u_inner))) * norm_ + 0) * norm_ end """ @@ -339,10 +347,11 @@ Should be used together with [`TreeMesh`](@ref). surface_flux_function, equations::CompressibleEulerEquations2D) # get the appropriate normal vector from the orientation + RealT = eltype(u_inner) if orientation == 1 - normal_direction = SVector(1, 0) + normal_direction = SVector(one(RealT), zero(RealT)) else # orientation == 2 - normal_direction = SVector(0, 1) + normal_direction = SVector(zero(RealT), one(RealT)) end # compute and return the flux using `boundary_condition_slip_wall` routine above @@ -380,7 +389,7 @@ end rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) if orientation == 1 f1 = rho_v1 f2 = rho_v1 * v1 + p @@ -434,21 +443,21 @@ The modification is in the energy flux to guarantee pressure equilibrium and was rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) # Average each factor of products in flux - rho_avg = 1 / 2 * (rho_ll + rho_rr) - v1_avg = 1 / 2 * (v1_ll + v1_rr) - v2_avg = 1 / 2 * (v2_ll + v2_rr) - p_avg = 1 / 2 * (p_ll + p_rr) - kin_avg = 1 / 2 * (v1_ll * v1_rr + v2_ll * v2_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + kin_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) # Calculate fluxes depending on orientation if orientation == 1 - pv1_avg = 1 / 2 * (p_ll * v1_rr + p_rr * v1_ll) + pv1_avg = 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) f1 = rho_avg * v1_avg f2 = f1 * v1_avg + p_avg f3 = f1 * v2_avg f4 = p_avg * v1_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv1_avg else - pv2_avg = 1 / 2 * (p_ll * v2_rr + p_rr * v2_ll) + pv2_avg = 0.5f0 * (p_ll * v2_rr + p_rr * v2_ll) f1 = rho_avg * v2_avg f2 = f1 * v1_avg f3 = f1 * v2_avg + p_avg @@ -467,12 +476,12 @@ end v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] # Average each factor of products in flux - rho_avg = 1 / 2 * (rho_ll + rho_rr) - v1_avg = 1 / 2 * (v1_ll + v1_rr) - v2_avg = 1 / 2 * (v2_ll + v2_rr) - v_dot_n_avg = 1 / 2 * (v_dot_n_ll + v_dot_n_rr) - p_avg = 1 / 2 * (p_ll + p_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v_dot_n_avg = 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) # Calculate fluxes depending on normal_direction f1 = rho_avg * v_dot_n_avg @@ -480,7 +489,7 @@ end f3 = f1 * v2_avg + p_avg * normal_direction[2] f4 = (f1 * velocity_square_avg + p_avg * v_dot_n_avg * equations.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)) return SVector(f1, f2, f3, f4) end @@ -504,11 +513,11 @@ Kinetic energy preserving two-point flux by rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) # Average each factor of products in flux - rho_avg = 1 / 2 * (rho_ll + rho_rr) - v1_avg = 1 / 2 * (v1_ll + v1_rr) - v2_avg = 1 / 2 * (v2_ll + v2_rr) - p_avg = 1 / 2 * (p_ll + p_rr) - e_avg = 1 / 2 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + e_avg = 0.5f0 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) # Calculate fluxes depending on orientation if orientation == 1 @@ -535,12 +544,12 @@ end rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) # Average each factor of products in flux - rho_avg = 0.5 * (rho_ll + rho_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) v_dot_n_avg = v1_avg * normal_direction[1] + v2_avg * normal_direction[2] - p_avg = 0.5 * (p_ll + p_rr) - e_avg = 0.5 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + e_avg = 0.5f0 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) # Calculate fluxes depending on normal_direction f1 = rho_avg * v_dot_n_avg @@ -565,19 +574,19 @@ Entropy conserving two-point flux by # Unpack left and right state rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) - beta_ll = 0.5 * rho_ll / p_ll - beta_rr = 0.5 * rho_rr / p_rr - specific_kin_ll = 0.5 * (v1_ll^2 + v2_ll^2) - specific_kin_rr = 0.5 * (v1_rr^2 + v2_rr^2) + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr + specific_kin_ll = 0.5f0 * (v1_ll^2 + v2_ll^2) + specific_kin_rr = 0.5f0 * (v1_rr^2 + v2_rr^2) # Compute the necessary mean values - rho_avg = 0.5 * (rho_ll + rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) rho_mean = ln_mean(rho_ll, rho_rr) beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_mean = 0.5 * rho_avg / beta_avg + beta_avg = 0.5f0 * (beta_ll + beta_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_mean = 0.5f0 * rho_avg / beta_avg velocity_square_avg = specific_kin_ll + specific_kin_rr # Calculate fluxes depending on orientation @@ -585,13 +594,15 @@ Entropy conserving two-point flux by f1 = rho_mean * v1_avg f2 = f1 * v1_avg + p_mean f3 = f1 * v2_avg - f4 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f4 = f1 * 0.5f0 * + (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + f2 * v1_avg + f3 * v2_avg else f1 = rho_mean * v2_avg f2 = f1 * v1_avg f3 = f1 * v2_avg + p_mean - f4 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f4 = f1 * 0.5f0 * + (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + f2 * v1_avg + f3 * v2_avg end @@ -605,26 +616,26 @@ end rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] - beta_ll = 0.5 * rho_ll / p_ll - beta_rr = 0.5 * rho_rr / p_rr - specific_kin_ll = 0.5 * (v1_ll^2 + v2_ll^2) - specific_kin_rr = 0.5 * (v1_rr^2 + v2_rr^2) + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr + specific_kin_ll = 0.5f0 * (v1_ll^2 + v2_ll^2) + specific_kin_rr = 0.5f0 * (v1_rr^2 + v2_rr^2) # Compute the necessary mean values - rho_avg = 0.5 * (rho_ll + rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) rho_mean = ln_mean(rho_ll, rho_rr) beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_mean = 0.5 * rho_avg / beta_avg + beta_avg = 0.5f0 * (beta_ll + beta_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_mean = 0.5f0 * rho_avg / beta_avg velocity_square_avg = specific_kin_ll + specific_kin_rr # Multiply with average of normal velocities - f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) f2 = f1 * v1_avg + p_mean * normal_direction[1] f3 = f1 * v2_avg + p_mean * normal_direction[2] - f4 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f4 = f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + f2 * v1_avg + f3 * v2_avg return SVector(f1, f2, f3, f4) @@ -658,10 +669,10 @@ See also # 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) - p_avg = 0.5 * (p_ll + p_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) # Calculate fluxes depending on orientation if orientation == 1 @@ -670,14 +681,14 @@ See also f3 = f1 * v2_avg f4 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.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 f1 = rho_mean * v2_avg f2 = f1 * v1_avg f3 = f1 * v2_avg + p_avg f4 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.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 return SVector(f1, f2, f3, f4) @@ -698,18 +709,18 @@ end # 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) - p_avg = 0.5 * (p_ll + p_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) # Calculate fluxes depending on normal_direction - f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) f2 = f1 * v1_avg + p_avg * normal_direction[1] f3 = f1 * v2_avg + p_avg * normal_direction[2] f4 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.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)) return SVector(f1, f2, f3, f4) end @@ -752,7 +763,7 @@ end rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) if orientation == 1 @@ -766,12 +777,12 @@ end alpha_p = 2 * (equations.gamma - 1) * lambda1_p + lambda2_p + lambda3_p - rho_2gamma = 0.5 * rho / equations.gamma + rho_2gamma = 0.5f0 * rho / equations.gamma f1p = rho_2gamma * alpha_p f2p = rho_2gamma * (alpha_p * v1 + a * (lambda2_p - lambda3_p)) f3p = rho_2gamma * alpha_p * v2 f4p = rho_2gamma * - (alpha_p * 0.5 * (v1^2 + v2^2) + a * v1 * (lambda2_p - lambda3_p) + (alpha_p * 0.5f0 * (v1^2 + v2^2) + a * v1 * (lambda2_p - lambda3_p) + a^2 * (lambda2_p + lambda3_p) * equations.inv_gamma_minus_one) else # orientation == 2 lambda1 = v2 @@ -784,12 +795,12 @@ end alpha_p = 2 * (equations.gamma - 1) * lambda1_p + lambda2_p + lambda3_p - rho_2gamma = 0.5 * rho / equations.gamma + rho_2gamma = 0.5f0 * rho / equations.gamma f1p = rho_2gamma * alpha_p f2p = rho_2gamma * alpha_p * v1 f3p = rho_2gamma * (alpha_p * v2 + a * (lambda2_p - lambda3_p)) f4p = rho_2gamma * - (alpha_p * 0.5 * (v1^2 + v2^2) + a * v2 * (lambda2_p - lambda3_p) + (alpha_p * 0.5f0 * (v1^2 + v2^2) + a * v2 * (lambda2_p - lambda3_p) + a^2 * (lambda2_p + lambda3_p) * equations.inv_gamma_minus_one) end return SVector(f1p, f2p, f3p, f4p) @@ -800,7 +811,7 @@ end rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) if orientation == 1 @@ -814,12 +825,12 @@ end alpha_m = 2 * (equations.gamma - 1) * lambda1_m + lambda2_m + lambda3_m - rho_2gamma = 0.5 * rho / equations.gamma + rho_2gamma = 0.5f0 * rho / equations.gamma f1m = rho_2gamma * alpha_m f2m = rho_2gamma * (alpha_m * v1 + a * (lambda2_m - lambda3_m)) f3m = rho_2gamma * alpha_m * v2 f4m = rho_2gamma * - (alpha_m * 0.5 * (v1^2 + v2^2) + a * v1 * (lambda2_m - lambda3_m) + (alpha_m * 0.5f0 * (v1^2 + v2^2) + a * v1 * (lambda2_m - lambda3_m) + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one) else # orientation == 2 lambda1 = v2 @@ -832,12 +843,12 @@ end alpha_m = 2 * (equations.gamma - 1) * lambda1_m + lambda2_m + lambda3_m - rho_2gamma = 0.5 * rho / equations.gamma + rho_2gamma = 0.5f0 * rho / equations.gamma f1m = rho_2gamma * alpha_m f2m = rho_2gamma * alpha_m * v1 f3m = rho_2gamma * (alpha_m * v2 + a * (lambda2_m - lambda3_m)) f4m = rho_2gamma * - (alpha_m * 0.5 * (v1^2 + v2^2) + a * v2 * (lambda2_m - lambda3_m) + (alpha_m * 0.5f0 * (v1^2 + v2^2) + a * v2 * (lambda2_m - lambda3_m) + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one) end return SVector(f1m, f2m, f3m, f4m) @@ -888,7 +899,7 @@ end rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho @@ -899,8 +910,8 @@ end lambda1_p = positive_part(lambda1) # Same as (lambda_i + abs(lambda_i)) / 2, but faster :) lambda2_p = positive_part(lambda2) - rhoa_2gamma = 0.5 * rho * a / equations.gamma - f1p = 0.5 * rho * (lambda1_p + lambda2_p) + rhoa_2gamma = 0.5f0 * rho * a / equations.gamma + f1p = 0.5f0 * rho * (lambda1_p + lambda2_p) f2p = f1p * v1 + rhoa_2gamma * (lambda1_p - lambda2_p) f3p = f1p * v2 f4p = f1p * H @@ -911,8 +922,8 @@ end lambda1_p = positive_part(lambda1) # Same as (lambda_i + abs(lambda_i)) / 2, but faster :) lambda2_p = positive_part(lambda2) - rhoa_2gamma = 0.5 * rho * a / equations.gamma - f1p = 0.5 * rho * (lambda1_p + lambda2_p) + rhoa_2gamma = 0.5f0 * rho * a / equations.gamma + f1p = 0.5f0 * rho * (lambda1_p + lambda2_p) f2p = f1p * v1 f3p = f1p * v2 + rhoa_2gamma * (lambda1_p - lambda2_p) f4p = f1p * H @@ -925,7 +936,7 @@ end rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho @@ -936,8 +947,8 @@ end lambda1_m = negative_part(lambda1) # Same as (lambda_i - abs(lambda_i)) / 2, but faster :) lambda2_m = negative_part(lambda2) - rhoa_2gamma = 0.5 * rho * a / equations.gamma - f1m = 0.5 * rho * (lambda1_m + lambda2_m) + rhoa_2gamma = 0.5f0 * rho * a / equations.gamma + f1m = 0.5f0 * rho * (lambda1_m + lambda2_m) f2m = f1m * v1 + rhoa_2gamma * (lambda1_m - lambda2_m) f3m = f1m * v2 f4m = f1m * H @@ -948,8 +959,8 @@ end lambda1_m = negative_part(lambda1) # Same as (lambda_i - abs(lambda_i)) / 2, but faster :) lambda2_m = negative_part(lambda2) - rhoa_2gamma = 0.5 * rho * a / equations.gamma - f1m = 0.5 * rho * (lambda1_m + lambda2_m) + rhoa_2gamma = 0.5f0 * rho * a / equations.gamma + f1m = 0.5f0 * rho * (lambda1_m + lambda2_m) f2m = f1m * v1 f3m = f1m * v2 + rhoa_2gamma * (lambda1_m - lambda2_m) f4m = f1m * H @@ -963,7 +974,7 @@ end rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho @@ -975,8 +986,8 @@ end lambda1_p = positive_part(lambda1) # Same as (lambda_i + abs(lambda_i)) / 2, but faster :) lambda2_p = positive_part(lambda2) - rhoa_2gamma = 0.5 * rho * a / equations.gamma - f1p = 0.5 * rho * (lambda1_p + lambda2_p) + rhoa_2gamma = 0.5f0 * rho * a / equations.gamma + f1p = 0.5f0 * rho * (lambda1_p + lambda2_p) f2p = f1p * v1 + rhoa_2gamma * normal_direction[1] * (lambda1_p - lambda2_p) f3p = f1p * v2 + rhoa_2gamma * normal_direction[2] * (lambda1_p - lambda2_p) f4p = f1p * H @@ -990,7 +1001,7 @@ end rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho @@ -1002,8 +1013,8 @@ end lambda1_m = negative_part(lambda1) # Same as (lambda_i - abs(lambda_i)) / 2, but faster :) lambda2_m = negative_part(lambda2) - rhoa_2gamma = 0.5 * rho * a / equations.gamma - f1m = 0.5 * rho * (lambda1_m + lambda2_m) + rhoa_2gamma = 0.5f0 * rho * a / equations.gamma + f1m = 0.5f0 * rho * (lambda1_m + lambda2_m) f2m = f1m * v1 + rhoa_2gamma * normal_direction[1] * (lambda1_m - lambda2_m) f3m = f1m * v2 + rhoa_2gamma * normal_direction[2] * (lambda1_m - lambda2_m) f4m = f1m * H @@ -1045,9 +1056,9 @@ end v_rr = v2_rr end - rho = 0.5 * (rho_ll + rho_rr) - p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) - v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) + rho = 0.5f0 * (rho_ll + rho_rr) + p = 0.5f0 * (p_ll + p_rr) - 0.5f0 * c * rho * (v_rr - v_ll) + v = 0.5f0 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) # We treat the energy term analogous to the potential temperature term in the paper by # Chen et al., i.e. we use p_ll and p_rr, and not p @@ -1083,9 +1094,9 @@ end # and then multiplying v by `norm_` again, but this version is slightly faster. norm_ = norm(normal_direction) - rho = 0.5 * (rho_ll + rho_rr) - p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_ - v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ + rho = 0.5f0 * (rho_ll + rho_rr) + p = 0.5f0 * (p_ll + p_rr) - 0.5f0 * c * rho * (v_rr - v_ll) / norm_ + v = 0.5f0 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ # We treat the energy term analogous to the potential temperature term in the paper by # Chen et al., i.e. we use p_ll and p_rr, and not p @@ -1155,24 +1166,24 @@ end rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho if orientation == 1 M = v1 / a - p_plus = 0.5 * (1 + equations.gamma * M) * p + p_plus = 0.5f0 * (1 + equations.gamma * M) * p - f1p = 0.25 * rho * a * (M + 1)^2 + f1p = 0.25f0 * rho * a * (M + 1)^2 f2p = f1p * v1 + p_plus f3p = f1p * v2 f4p = f1p * H else # orientation == 2 M = v2 / a - p_plus = 0.5 * (1 + equations.gamma * M) * p + p_plus = 0.5f0 * (1 + equations.gamma * M) * p - f1p = 0.25 * rho * a * (M + 1)^2 + f1p = 0.25f0 * rho * a * (M + 1)^2 f2p = f1p * v1 f3p = f1p * v2 + p_plus f4p = f1p * H @@ -1185,24 +1196,24 @@ end rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho if orientation == 1 M = v1 / a - p_minus = 0.5 * (1 - equations.gamma * M) * p + p_minus = 0.5f0 * (1 - equations.gamma * M) * p - f1m = -0.25 * rho * a * (M - 1)^2 + f1m = -0.25f0 * rho * a * (M - 1)^2 f2m = f1m * v1 + p_minus f3m = f1m * v2 f4m = f1m * H else # orientation == 2 M = v2 / a - p_minus = 0.5 * (1 - equations.gamma * M) * p + p_minus = 0.5f0 * (1 - equations.gamma * M) * p - f1m = -0.25 * rho * a * (M - 1)^2 + f1m = -0.25f0 * rho * a * (M - 1)^2 f2m = f1m * v1 f3m = f1m * v2 + p_minus f4m = f1m * H @@ -1216,16 +1227,16 @@ end rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho v_n = normal_direction[1] * v1 + normal_direction[2] * v2 M = v_n / a - p_plus = 0.5 * (1 + equations.gamma * M) * p + p_plus = 0.5f0 * (1 + equations.gamma * M) * p - f1p = 0.25 * rho * a * (M + 1)^2 + f1p = 0.25f0 * rho * a * (M + 1)^2 f2p = f1p * v1 + normal_direction[1] * p_plus f3p = f1p * v2 + normal_direction[2] * p_plus f4p = f1p * H @@ -1239,16 +1250,16 @@ end rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho v_n = normal_direction[1] * v1 + normal_direction[2] * v2 M = v_n / a - p_minus = 0.5 * (1 - equations.gamma * M) * p + p_minus = 0.5f0 * (1 - equations.gamma * M) * p - f1m = -0.25 * rho * a * (M - 1)^2 + f1m = -0.25f0 * rho * a * (M - 1)^2 f2m = f1m * v1 + normal_direction[1] * p_minus f3m = f1m * v2 + normal_direction[2] * p_minus f4m = f1m * H @@ -1289,24 +1300,24 @@ end rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho - lambda = 0.5 * (sqrt(v1^2 + v2^2) + a) + lambda = 0.5f0 * (sqrt(v1^2 + v2^2) + a) if orientation == 1 #lambda = 0.5 * (abs(v1) + a) - f1p = 0.5 * rho * v1 + lambda * u[1] - f2p = 0.5 * rho * v1 * v1 + 0.5 * p + lambda * u[2] - f3p = 0.5 * rho * v1 * v2 + lambda * u[3] - f4p = 0.5 * rho * v1 * H + lambda * u[4] + f1p = 0.5f0 * rho * v1 + lambda * u[1] + f2p = 0.5f0 * rho * v1 * v1 + 0.5f0 * p + lambda * u[2] + f3p = 0.5f0 * rho * v1 * v2 + lambda * u[3] + f4p = 0.5f0 * rho * v1 * H + lambda * u[4] else # orientation == 2 #lambda = 0.5 * (abs(v2) + a) - f1p = 0.5 * rho * v2 + lambda * u[1] - f2p = 0.5 * rho * v2 * v1 + lambda * u[2] - f3p = 0.5 * rho * v2 * v2 + 0.5 * p + lambda * u[3] - f4p = 0.5 * rho * v2 * H + lambda * u[4] + f1p = 0.5f0 * rho * v2 + lambda * u[1] + f2p = 0.5f0 * rho * v2 * v1 + lambda * u[2] + f3p = 0.5f0 * rho * v2 * v2 + 0.5f0 * p + lambda * u[3] + f4p = 0.5f0 * rho * v2 * H + lambda * u[4] end return SVector(f1p, f2p, f3p, f4p) end @@ -1316,24 +1327,24 @@ end rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho - lambda = 0.5 * (sqrt(v1^2 + v2^2) + a) + lambda = 0.5f0 * (sqrt(v1^2 + v2^2) + a) if orientation == 1 #lambda = 0.5 * (abs(v1) + a) - f1m = 0.5 * rho * v1 - lambda * u[1] - f2m = 0.5 * rho * v1 * v1 + 0.5 * p - lambda * u[2] - f3m = 0.5 * rho * v1 * v2 - lambda * u[3] - f4m = 0.5 * rho * v1 * H - lambda * u[4] + f1m = 0.5f0 * rho * v1 - lambda * u[1] + f2m = 0.5f0 * rho * v1 * v1 + 0.5f0 * p - lambda * u[2] + f3m = 0.5f0 * rho * v1 * v2 - lambda * u[3] + f4m = 0.5f0 * rho * v1 * H - lambda * u[4] else # orientation == 2 #lambda = 0.5 * (abs(v2) + a) - f1m = 0.5 * rho * v2 - lambda * u[1] - f2m = 0.5 * rho * v2 * v1 - lambda * u[2] - f3m = 0.5 * rho * v2 * v2 + 0.5 * p - lambda * u[3] - f4m = 0.5 * rho * v2 * H - lambda * u[4] + f1m = 0.5f0 * rho * v2 - lambda * u[1] + f2m = 0.5f0 * rho * v2 * v1 - lambda * u[2] + f3m = 0.5f0 * rho * v2 * v2 + 0.5f0 * p - lambda * u[3] + f4m = 0.5f0 * rho * v2 * H - lambda * u[4] end return SVector(f1m, f2m, f3m, f4m) end @@ -1346,15 +1357,15 @@ end a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho - lambda = 0.5 * (sqrt(v1^2 + v2^2) + a) + lambda = 0.5f0 * (sqrt(v1^2 + v2^2) + a) v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] rho_v_normal = rho * v_normal - f1p = 0.5 * rho_v_normal + lambda * u[1] - f2p = 0.5 * rho_v_normal * v1 + 0.5 * p * normal_direction[1] + lambda * u[2] - f3p = 0.5 * rho_v_normal * v2 + 0.5 * p * normal_direction[2] + lambda * u[3] - f4p = 0.5 * rho_v_normal * H + lambda * u[4] + f1p = 0.5f0 * rho_v_normal + lambda * u[1] + f2p = 0.5f0 * rho_v_normal * v1 + 0.5f0 * p * normal_direction[1] + lambda * u[2] + f3p = 0.5f0 * rho_v_normal * v2 + 0.5f0 * p * normal_direction[2] + lambda * u[3] + f4p = 0.5f0 * rho_v_normal * H + lambda * u[4] return SVector(f1p, f2p, f3p, f4p) end @@ -1367,15 +1378,15 @@ end a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho - lambda = 0.5 * (sqrt(v1^2 + v2^2) + a) + lambda = 0.5f0 * (sqrt(v1^2 + v2^2) + a) v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] rho_v_normal = rho * v_normal - f1m = 0.5 * rho_v_normal - lambda * u[1] - f2m = 0.5 * rho_v_normal * v1 + 0.5 * p * normal_direction[1] - lambda * u[2] - f3m = 0.5 * rho_v_normal * v2 + 0.5 * p * normal_direction[2] - lambda * u[3] - f4m = 0.5 * rho_v_normal * H - lambda * u[4] + f1m = 0.5f0 * rho_v_normal - lambda * u[1] + f2m = 0.5f0 * rho_v_normal * v1 + 0.5f0 * p * normal_direction[1] - lambda * u[2] + f3m = 0.5f0 * rho_v_normal * v2 + 0.5f0 * p * normal_direction[2] - lambda * u[3] + f4m = 0.5f0 * rho_v_normal * H - lambda * u[4] return SVector(f1m, f2m, f3m, f4m) end @@ -1555,13 +1566,13 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, v1_ll = rho_v1_ll / rho_ll v2_ll = rho_v2_ll / rho_ll e_ll = rho_e_ll / rho_ll - p_ll = (equations.gamma - 1) * (rho_e_ll - 1 / 2 * rho_ll * (v1_ll^2 + v2_ll^2)) + p_ll = (equations.gamma - 1) * (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2 + v2_ll^2)) c_ll = sqrt(equations.gamma * p_ll / rho_ll) v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr e_rr = rho_e_rr / rho_rr - p_rr = (equations.gamma - 1) * (rho_e_rr - 1 / 2 * rho_rr * (v1_rr^2 + v2_rr^2)) + p_rr = (equations.gamma - 1) * (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2 + v2_rr^2)) c_rr = sqrt(equations.gamma * p_rr / rho_rr) # Obtain left and right fluxes @@ -1586,18 +1597,18 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, H_ll = (rho_e_ll + p_ll) / rho_ll H_rr = (rho_e_rr + p_rr) / rho_rr H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5f0 * vel_roe_mag)) Ssl = min(vel_L - c_ll, vel_roe - c_roe) Ssr = max(vel_R + c_rr, vel_roe + c_roe) sMu_L = Ssl - vel_L sMu_R = Ssr - vel_R - if Ssl >= 0.0 + if Ssl >= 0 f1 = f_ll[1] f2 = f_ll[2] f3 = f_ll[3] f4 = f_ll[4] - elseif Ssr <= 0.0 + elseif Ssr <= 0 f1 = f_rr[1] f2 = f_rr[2] f3 = f_rr[3] @@ -1605,7 +1616,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, else SStar = (p_rr - p_ll + rho_ll * vel_L * sMu_L - rho_rr * vel_R * sMu_R) / (rho_ll * sMu_L - rho_rr * sMu_R) - if Ssl <= 0.0 <= SStar + if Ssl <= 0 <= SStar densStar = rho_ll * sMu_L / (Ssl - SStar) enerStar = e_ll + (SStar - vel_L) * (SStar + p_ll / (rho_ll * sMu_L)) UStar1 = densStar @@ -1679,19 +1690,19 @@ function flux_hllc(u_ll, u_rr, normal_direction::AbstractVector, H_rr = (u_rr[4] + p_rr) / rho_rr H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) * norm_ + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5f0 * vel_roe_mag)) * norm_ Ssl = min(v_dot_n_ll - c_ll, vel_roe - c_roe) Ssr = max(v_dot_n_rr + c_rr, vel_roe + c_roe) sMu_L = Ssl - v_dot_n_ll sMu_R = Ssr - v_dot_n_rr - if Ssl >= 0.0 + if Ssl >= 0 f1 = f_ll[1] f2 = f_ll[2] f3 = f_ll[3] f4 = f_ll[4] - elseif Ssr <= 0.0 + elseif Ssr <= 0 f1 = f_rr[1] f2 = f_rr[2] f3 = f_rr[3] @@ -1699,7 +1710,7 @@ function flux_hllc(u_ll, u_rr, normal_direction::AbstractVector, else SStar = (rho_ll * v_dot_n_ll * sMu_L - rho_rr * v_dot_n_rr * sMu_R + (p_rr - p_ll) * norm_sq) / (rho_ll * sMu_L - rho_rr * sMu_R) - if Ssl <= 0.0 <= SStar + if Ssl <= 0 <= SStar densStar = rho_ll * sMu_L / (Ssl - SStar) enerStar = e_ll + (SStar - v_dot_n_ll) * @@ -1773,19 +1784,19 @@ of the numerical flux. v_roe_mag = v1_roe^2 + v2_roe^2 H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) * inv_sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag)) + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5f0 * v_roe_mag)) # Compute convenience constant for positivity preservation, see # https://doi.org/10.1016/0021-9991(91)90211-3 - beta = sqrt(0.5 * (equations.gamma - 1) / equations.gamma) + beta = sqrt(0.5f0 * (equations.gamma - 1) / equations.gamma) # Estimate the edges of the Riemann fan (with positivity conservation) if orientation == 1 # x-direction - SsL = min(v1_roe - c_roe, v1_ll - beta * c_ll, zero(v1_roe)) - SsR = max(v1_roe + c_roe, v1_rr + beta * c_rr, zero(v1_roe)) + SsL = min(v1_roe - c_roe, v1_ll - beta * c_ll, 0) + SsR = max(v1_roe + c_roe, v1_rr + beta * c_rr, 0) elseif orientation == 2 # y-direction - SsL = min(v2_roe - c_roe, v2_ll - beta * c_ll, zero(v2_roe)) - SsR = max(v2_roe + c_roe, v2_rr + beta * c_rr, zero(v2_roe)) + SsL = min(v2_roe - c_roe, v2_ll - beta * c_ll, 0) + SsR = max(v2_roe + c_roe, v2_rr + beta * c_rr, 0) end return SsL, SsR @@ -1836,15 +1847,15 @@ of the numerical flux. v_roe_mag = v1_roe^2 + v2_roe^2 H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) * inv_sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag)) * norm_ + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5f0 * v_roe_mag)) * norm_ # Compute convenience constant for positivity preservation, see # https://doi.org/10.1016/0021-9991(91)90211-3 - beta = sqrt(0.5 * (equations.gamma - 1) / equations.gamma) + beta = sqrt(0.5f0 * (equations.gamma - 1) / equations.gamma) # Estimate the edges of the Riemann fan (with positivity conservation) - SsL = min(v_roe - c_roe, v_dot_n_ll - beta * c_ll, zero(v_roe)) - SsR = max(v_roe + c_roe, v_dot_n_rr + beta * c_rr, zero(v_roe)) + SsL = min(v_roe - c_roe, v_dot_n_ll - beta * c_ll, 0) + SsR = max(v_roe + c_roe, v_dot_n_rr + beta * c_rr, 0) return SsL, SsR end @@ -1862,7 +1873,7 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) return SVector(rho, v1, v2, p) end @@ -1874,11 +1885,12 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v_square = v1^2 + v2^2 - p = (equations.gamma - 1) * (rho_e - 0.5 * rho * v_square) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho * v_square) s = log(p) - equations.gamma * log(rho) rho_p = rho / p - w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - 0.5 * rho_p * v_square + w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - + 0.5f0 * rho_p * v_square w2 = rho_p * v1 w3 = rho_p * v2 w4 = -rho_p @@ -1895,11 +1907,11 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v_square = v1^2 + v2^2 - inv_rho_gammap1 = (1 / rho)^(equations.gamma + 1.0) + inv_rho_gammap1 = (1 / rho)^(equations.gamma + 1) # The derivative vector for the modified specific entropy of Guermond et al. w1 = inv_rho_gammap1 * - (0.5 * rho * (equations.gamma + 1.0) * v_square - equations.gamma * rho_e) + (0.5f0 * rho * (equations.gamma + 1) * v_square - equations.gamma * rho_e) w2 = -rho_v1 * inv_rho_gammap1 w3 = -rho_v2 * inv_rho_gammap1 w4 = (1 / rho)^equations.gamma @@ -1936,7 +1948,7 @@ end rho, v1, v2, p = prim rho_v1 = rho * v1 rho_v2 = rho * v2 - rho_e = p * equations.inv_gamma_minus_one + 0.5 * (rho_v1 * v1 + rho_v2 * v2) + rho_e = p * equations.inv_gamma_minus_one + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2) return SVector(rho, rho_v1, rho_v2, rho_e) end @@ -1947,7 +1959,7 @@ end @inline function pressure(u, equations::CompressibleEulerEquations2D) rho, rho_v1, rho_v2, rho_e = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2) / rho) return p end @@ -1960,12 +1972,12 @@ end v2 = rho_v2 / rho v_square = v1^2 + v2^2 - return (equations.gamma - 1.0) * SVector(0.5 * v_square, -v1, -v2, 1.0) + return (equations.gamma - 1) * SVector(0.5f0 * v_square, -v1, -v2, 1) end @inline function density_pressure(u, equations::CompressibleEulerEquations2D) rho, rho_v1, rho_v2, rho_e = u - rho_times_p = (equations.gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2)) + rho_times_p = (equations.gamma - 1) * (rho * rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2)) return rho_times_p end @@ -1995,7 +2007,7 @@ end # Calculate thermodynamic entropy for a conservative state `cons` @inline function entropy_thermodynamic(cons, equations::CompressibleEulerEquations2D) # Pressure - p = (equations.gamma - 1) * (cons[4] - 1 / 2 * (cons[2]^2 + cons[3]^2) / cons[1]) + p = (equations.gamma - 1) * (cons[4] - 0.5f0 * (cons[2]^2 + cons[3]^2) / cons[1]) # Thermodynamic entropy s = log(p) - equations.gamma * log(cons[1]) @@ -2034,7 +2046,7 @@ Note: This is *not* the "conventional" specific entropy ``s = ln(p / \rho^\gamma rho, rho_v1, rho_v2, rho_e = u # Modified specific entropy from Guermond et al. (2019) - s = (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho) * (1 / rho)^equations.gamma + s = (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2) / rho) * (1 / rho)^equations.gamma return s end @@ -2067,7 +2079,7 @@ end # State validation for Newton-bisection method of subcell IDP limiting @inline function Base.isvalid(u, equations::CompressibleEulerEquations2D) p = pressure(u, equations) - if u[1] <= 0.0 || p <= 0.0 + if u[1] <= 0 || p <= 0 return false end return true diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index f156aa29689..4f4d4553a8f 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -67,11 +67,12 @@ end A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equations::CompressibleEulerEquations3D) - rho = 1.0 - rho_v1 = 0.1 - rho_v2 = -0.2 - rho_v3 = 0.7 - rho_e = 10.0 + RealT = eltype(x) + rho = 1 + rho_v1 = convert(RealT, 0.1) + rho_v2 = convert(RealT, -0.2) + rho_v3 = convert(RealT, 0.7) + rho_e = 10 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e) end @@ -83,11 +84,12 @@ A smooth initial condition used for convergence tests in combination with """ function initial_condition_convergence_test(x, t, equations::CompressibleEulerEquations3D) + 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] + x[2] + x[3] - t)) rho = ini @@ -108,11 +110,12 @@ Source terms used for convergence tests in combination with @inline function source_terms_convergence_test(u, x, t, equations::CompressibleEulerEquations3D) # Same settings as in `initial_condition` + 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 γ = equations.gamma x1, x2, x3 = x @@ -121,7 +124,7 @@ Source terms used for convergence tests in combination with rho_x = ω * A * co # Note that d/dt rho = -d/dx rho = -d/dy rho = - d/dz rho. - tmp = (2 * rho - 1.5) * (γ - 1) + tmp = (2 * rho - 1.5f0) * (γ - 1) du1 = 2 * rho_x du2 = rho_x * (2 + tmp) @@ -144,20 +147,21 @@ function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations3D) # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # Set up spherical coordinates + RealT = eltype(x) inicenter = (0, 0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] z_norm = x[3] - inicenter[3] r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) phi = atan(y_norm, x_norm) - theta = iszero(r) ? 0.0 : acos(z_norm / r) + theta = iszero(r) ? zero(RealT) : acos(z_norm / r) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi) * sin(theta) - v2 = r > 0.5 ? 0.0 : 0.1882 * sin(phi) * sin(theta) - v3 = r > 0.5 ? 0.0 : 0.1882 * cos(theta) - p = r > 0.5 ? 1.0 : 1.245 + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi) * sin(theta) + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin(phi) * sin(theta) + v3 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(theta) + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) return prim2cons(SVector(rho, v1, v2, v3, p), equations) end @@ -175,19 +179,20 @@ or [`source_terms_eoc_test_euler`](@ref). function initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::CompressibleEulerEquations3D) # OBS! this assumes that γ = 2 other manufactured source terms are incorrect - if equations.gamma != 2.0 + if equations.gamma != 2 error("adiabatic constant must be 2 for the coupling convergence test") end - c = 2.0 - A = 0.1 - ini = c + A * sin(pi * (x[1] + x[2] + x[3] - t)) - G = 1.0 # gravitational constant + RealT = eltype(x) + c = 2 + A = convert(RealT, 0.1) + ini = c + A * sin(convert(RealT, pi) * (x[1] + x[2] + x[3] - t)) + G = 1 # gravitational constant rho = ini - v1 = 1.0 - v2 = 1.0 - v3 = 1.0 - p = ini^2 * G * 2 / (3 * pi) # "3" is the number of spatial dimensions + v1 = 1 + v2 = 1 + v3 = 1 + p = ini^2 * G * 2 / (3 * convert(RealT, pi)) # "3" is the number of spatial dimensions return prim2cons(SVector(rho, v1, v2, v3, p), equations) end @@ -204,15 +209,16 @@ in combination with [`initial_condition_eoc_test_coupled_euler_gravity`](@ref). @inline function source_terms_eoc_test_coupled_euler_gravity(u, x, t, equations::CompressibleEulerEquations3D) # Same settings as in `initial_condition_eoc_test_coupled_euler_gravity` - c = 2.0 - A = 0.1 - G = 1.0 # gravitational constant, must match coupling solver - C_grav = -4 * G / (3 * pi) # "3" is the number of spatial dimensions # 2D: -2.0*G/pi + RealT = eltype(u) + c = 2 + A = convert(RealT, 0.1) + G = 1 # gravitational constant, must match coupling solver + C_grav = -4 * G / (3 * convert(RealT, pi)) # "3" is the number of spatial dimensions # 2D: -2.0*G/pi x1, x2, x3 = x # TODO: sincospi - si, co = sincos(pi * (x1 + x2 + x3 - t)) - rhox = A * pi * co + si, co = sincos(convert(RealT, pi) * (x1 + x2 + x3 - t)) + rhox = A * convert(RealT, pi) * co rho = c + A * si # In "2 * rhox", the "2" is "number of spatial dimensions minus one" @@ -220,7 +226,7 @@ in combination with [`initial_condition_eoc_test_coupled_euler_gravity`](@ref). du2 = 2 * rhox du3 = 2 * rhox du4 = 2 * rhox - du5 = 2 * rhox * (3 / 2 - C_grav * rho) # "3" in "3/2" is the number of spatial dimensions + du5 = 2 * rhox * (1.5f0 - C_grav * rho) # "3" in "3/2" is the number of spatial dimensions return SVector(du1, du2, du3, du4, du5) end @@ -241,15 +247,16 @@ in combination with [`initial_condition_eoc_test_coupled_euler_gravity`](@ref). """ function source_terms_eoc_test_euler(u, x, t, equations::CompressibleEulerEquations3D) # Same settings as in `initial_condition_eoc_test_coupled_euler_gravity` - c = 2.0 - A = 0.1 - G = 1.0 - C_grav = -4 * G / (3 * pi) # "3" is the number of spatial dimensions + RealT = eltype(u) + c = 2 + A = convert(RealT, 0.1) + G = 1 + C_grav = -4 * G / (3 * convert(RealT, pi)) # "3" is the number of spatial dimensions x1, x2, x3 = x # TODO: sincospi - si, co = sincos(pi * (x1 + x2 + x3 - t)) - rhox = A * pi * co + si, co = sincos(convert(RealT, pi) * (x1 + x2 + x3 - t)) + rhox = A * convert(RealT, pi) * co rho = c + A * si du1 = rhox * 2 @@ -309,26 +316,26 @@ Details about the 1D pressure Riemann solution can be found in Section 6.3.3 of # Eleuterio F. Toro (2009) # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) - if v_normal <= 0.0 + if v_normal <= 0 sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed p_star = p_local * - (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * - equations.gamma * - equations.inv_gamma_minus_one) - else # v_normal > 0.0 + (1 + 0.5f0 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * + equations.gamma * + equations.inv_gamma_minus_one) + else # v_normal > 0 A = 2 / ((equations.gamma + 1) * rho_local) B = p_local * (equations.gamma - 1) / (equations.gamma + 1) p_star = p_local + - 0.5 * v_normal / A * + 0.5f0 * v_normal / A * (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) end # For the slip wall we directly set the flux as the normal velocity is zero - return SVector(zero(eltype(u_inner)), + return SVector(0, p_star * normal[1], p_star * normal[2], p_star * normal[3], - zero(eltype(u_inner))) * norm_ + 0) * norm_ end """ @@ -342,12 +349,13 @@ Should be used together with [`TreeMesh`](@ref). surface_flux_function, equations::CompressibleEulerEquations3D) # get the appropriate normal vector from the orientation + RealT = eltype(u_inner) if orientation == 1 - normal_direction = SVector(1.0, 0.0, 0.0) + normal_direction = SVector(one(RealT), zero(RealT), zero(RealT)) elseif orientation == 2 - normal_direction = SVector(0.0, 1.0, 0.0) + normal_direction = SVector(zero(RealT), one(RealT), zero(RealT)) else # orientation == 3 - normal_direction = SVector(0.0, 0.0, 1.0) + normal_direction = SVector(zero(RealT), zero(RealT), one(RealT)) end # compute and return the flux using `boundary_condition_slip_wall` routine above @@ -387,7 +395,7 @@ end v2 = rho_v2 / rho v3 = rho_v3 / rho p = (equations.gamma - 1) * - (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3)) + (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3)) if orientation == 1 f1 = rho_v1 f2 = rho_v1 * v1 + p @@ -410,17 +418,18 @@ end return SVector(f1, f2, f3, f4, f5) end -@inline function flux(u, normal::AbstractVector, +@inline function flux(u, normal_direction::AbstractVector, equations::CompressibleEulerEquations3D) rho_e = last(u) rho, v1, v2, v3, p = cons2prim(u, equations) - v_normal = v1 * normal[1] + v2 * normal[2] + v3 * normal[3] + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + + v3 * normal_direction[3] rho_v_normal = rho * v_normal f1 = rho_v_normal - f2 = rho_v_normal * v1 + p * normal[1] - f3 = rho_v_normal * v2 + p * normal[2] - f4 = rho_v_normal * v3 + p * normal[3] + f2 = rho_v_normal * v1 + p * normal_direction[1] + f3 = rho_v_normal * v2 + p * normal_direction[2] + f4 = rho_v_normal * v3 + p * normal_direction[3] f5 = (rho_e + p) * v_normal return SVector(f1, f2, f3, f4, f5) end @@ -448,30 +457,30 @@ The modification is in the energy flux to guarantee pressure equilibrium and was rho_rr, v1_rr, v2_rr, v3_rr, p_rr = cons2prim(u_rr, equations) # Average each factor of products in flux - rho_avg = 1 / 2 * (rho_ll + rho_rr) - v1_avg = 1 / 2 * (v1_ll + v1_rr) - v2_avg = 1 / 2 * (v2_ll + v2_rr) - v3_avg = 1 / 2 * (v3_ll + v3_rr) - p_avg = 1 / 2 * (p_ll + p_rr) - kin_avg = 1 / 2 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + rho_avg = 0.5f0 * (rho_ll + rho_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) + kin_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) # Calculate fluxes depending on orientation if orientation == 1 - pv1_avg = 1 / 2 * (p_ll * v1_rr + p_rr * v1_ll) + pv1_avg = 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) f1 = rho_avg * v1_avg f2 = f1 * v1_avg + p_avg f3 = f1 * v2_avg f4 = f1 * v3_avg f5 = p_avg * v1_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv1_avg elseif orientation == 2 - pv2_avg = 1 / 2 * (p_ll * v2_rr + p_rr * v2_ll) + pv2_avg = 0.5f0 * (p_ll * v2_rr + p_rr * v2_ll) f1 = rho_avg * v2_avg f2 = f1 * v1_avg f3 = f1 * v2_avg + p_avg f4 = f1 * v3_avg f5 = p_avg * v2_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv2_avg else - pv3_avg = 1 / 2 * (p_ll * v3_rr + p_rr * v3_ll) + pv3_avg = 0.5f0 * (p_ll * v3_rr + p_rr * v3_ll) f1 = rho_avg * v3_avg f2 = f1 * v1_avg f3 = f1 * v2_avg @@ -493,13 +502,13 @@ end v3_rr * normal_direction[3] # Average each factor of products in flux - rho_avg = 1 / 2 * (rho_ll + rho_rr) - v1_avg = 1 / 2 * (v1_ll + v1_rr) - v2_avg = 1 / 2 * (v2_ll + v2_rr) - v3_avg = 1 / 2 * (v3_ll + v3_rr) - v_dot_n_avg = 1 / 2 * (v_dot_n_ll + v_dot_n_rr) - p_avg = 1 / 2 * (p_ll + p_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + rho_avg = 0.5f0 * (rho_ll + rho_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_dot_n_avg = 0.5f0 * (v_dot_n_ll + v_dot_n_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) # Calculate fluxes depending on normal_direction f1 = rho_avg * v_dot_n_avg @@ -508,7 +517,7 @@ end f4 = f1 * v3_avg + p_avg * normal_direction[3] f5 = (f1 * velocity_square_avg + p_avg * v_dot_n_avg * equations.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)) return SVector(f1, f2, f3, f4, f5) end @@ -532,12 +541,12 @@ Kinetic energy preserving two-point flux by rho_rr, v1_rr, v2_rr, v3_rr, p_rr = cons2prim(u_rr, equations) # Average each factor of products in flux - rho_avg = 0.5 * (rho_ll + rho_rr) - 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) - e_avg = 0.5 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_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) + e_avg = 0.5f0 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) # Calculate fluxes depending on orientation if orientation == 1 @@ -579,17 +588,17 @@ end v3_rr = rho_v3_rr / rho_rr # Average each factor of products in flux - rho_avg = 0.5 * (rho_ll + rho_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) + rho_avg = 0.5f0 * (rho_ll + rho_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_dot_n_avg = v1_avg * normal_direction[1] + v2_avg * normal_direction[2] + v3_avg * normal_direction[3] - p_avg = 0.5 * ((equations.gamma - 1) * - (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2 + v3_ll^2)) + + p_avg = 0.5f0 * ((equations.gamma - 1) * + (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2 + v2_ll^2 + v3_ll^2)) + (equations.gamma - 1) * - (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2 + v3_rr^2))) - e_avg = 0.5 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2 + v2_rr^2 + v3_rr^2))) + e_avg = 0.5f0 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) # Calculate fluxes depending on normal_direction f1 = rho_avg * v_dot_n_avg @@ -616,20 +625,20 @@ Entropy conserving two-point flux by rho_ll, v1_ll, v2_ll, v3_ll, p_ll = cons2prim(u_ll, equations) rho_rr, v1_rr, v2_rr, v3_rr, p_rr = cons2prim(u_rr, equations) - beta_ll = 0.5 * rho_ll / p_ll - beta_rr = 0.5 * rho_rr / p_rr - specific_kin_ll = 0.5 * (v1_ll^2 + v2_ll^2 + v3_ll^2) - specific_kin_rr = 0.5 * (v1_rr^2 + v2_rr^2 + v3_rr^2) + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr + specific_kin_ll = 0.5f0 * (v1_ll^2 + v2_ll^2 + v3_ll^2) + specific_kin_rr = 0.5f0 * (v1_rr^2 + v2_rr^2 + v3_rr^2) # Compute the necessary mean values - rho_avg = 0.5 * (rho_ll + rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) rho_mean = ln_mean(rho_ll, rho_rr) beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - 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_mean = 0.5 * rho_avg / beta_avg + beta_avg = 0.5f0 * (beta_ll + beta_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_mean = 0.5f0 * rho_avg / beta_avg velocity_square_avg = specific_kin_ll + specific_kin_rr # Calculate fluxes depending on orientation @@ -638,21 +647,24 @@ Entropy conserving two-point flux by f2 = f1 * v1_avg + p_mean f3 = f1 * v2_avg f4 = f1 * v3_avg - f5 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f5 = f1 * 0.5f0 * + (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg elseif orientation == 2 f1 = rho_mean * v2_avg f2 = f1 * v1_avg f3 = f1 * v2_avg + p_mean f4 = f1 * v3_avg - f5 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f5 = f1 * 0.5f0 * + (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg else f1 = rho_mean * v3_avg f2 = f1 * v1_avg f3 = f1 * v2_avg f4 = f1 * v3_avg + p_mean - f5 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f5 = f1 * 0.5f0 * + (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg end @@ -670,28 +682,28 @@ end v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + v3_rr * normal_direction[3] - beta_ll = 0.5 * rho_ll / p_ll - beta_rr = 0.5 * rho_rr / p_rr - specific_kin_ll = 0.5 * (v1_ll^2 + v2_ll^2 + v3_ll^2) - specific_kin_rr = 0.5 * (v1_rr^2 + v2_rr^2 + v3_rr^2) + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr + specific_kin_ll = 0.5f0 * (v1_ll^2 + v2_ll^2 + v3_ll^2) + specific_kin_rr = 0.5f0 * (v1_rr^2 + v2_rr^2 + v3_rr^2) # Compute the necessary mean values - rho_avg = 0.5 * (rho_ll + rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) rho_mean = ln_mean(rho_ll, rho_rr) beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - 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_mean = 0.5 * rho_avg / beta_avg + beta_avg = 0.5f0 * (beta_ll + beta_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_mean = 0.5f0 * rho_avg / beta_avg velocity_square_avg = specific_kin_ll + specific_kin_rr # Multiply with average of normal velocities - f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) f2 = f1 * v1_avg + p_mean * normal_direction[1] f3 = f1 * v2_avg + p_mean * normal_direction[2] f4 = f1 * v3_avg + p_mean * normal_direction[3] - f5 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f5 = f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg return SVector(f1, f2, f3, f4, f5) @@ -725,11 +737,11 @@ See also # 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) + 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) # Calculate fluxes depending on orientation if orientation == 1 @@ -739,7 +751,7 @@ See also f4 = f1 * v3_avg f5 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (p_ll * v1_rr + p_rr * v1_ll) + 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) elseif orientation == 2 f1 = rho_mean * v2_avg f2 = f1 * v1_avg @@ -747,7 +759,7 @@ See also f4 = f1 * v3_avg f5 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (p_ll * v2_rr + p_rr * v2_ll) + 0.5f0 * (p_ll * v2_rr + p_rr * v2_ll) else # orientation == 3 f1 = rho_mean * v3_avg f2 = f1 * v1_avg @@ -755,7 +767,7 @@ See also f4 = f1 * v3_avg + p_avg f5 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (p_ll * v3_rr + p_rr * v3_ll) + 0.5f0 * (p_ll * v3_rr + p_rr * v3_ll) end return SVector(f1, f2, f3, f4, f5) @@ -778,20 +790,20 @@ end # 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) + 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) # Calculate fluxes depending on normal_direction - f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) f2 = f1 * v1_avg + p_avg * normal_direction[1] f3 = f1 * v2_avg + p_avg * normal_direction[2] f4 = f1 * v3_avg + p_avg * normal_direction[3] f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.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)) return SVector(f1, f2, f3, f4, f5) end @@ -834,7 +846,7 @@ end v2 = rho_v2 / rho v3 = rho_v3 / rho p = (equations.gamma - 1) * - (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3)) + (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3)) a = sqrt(equations.gamma * p / rho) if orientation == 1 @@ -848,13 +860,15 @@ end alpha_p = 2 * (equations.gamma - 1) * lambda1_p + lambda2_p + lambda3_p - rho_2gamma = 0.5 * rho / equations.gamma + rho_2gamma = 0.5f0 * rho / equations.gamma f1p = rho_2gamma * alpha_p f2p = rho_2gamma * (alpha_p * v1 + a * (lambda2_p - lambda3_p)) f3p = rho_2gamma * alpha_p * v2 f4p = rho_2gamma * alpha_p * v3 f5p = rho_2gamma * - (alpha_p * 0.5 * (v1^2 + v2^2 + v3^2) + a * v1 * (lambda2_p - lambda3_p) + (alpha_p * 0.5f0 * (v1^2 + v2^2 + v3^2) + + a * v1 * + (lambda2_p - lambda3_p) + a^2 * (lambda2_p + lambda3_p) * equations.inv_gamma_minus_one) elseif orientation == 2 lambda1 = v2 @@ -867,13 +881,15 @@ end alpha_p = 2 * (equations.gamma - 1) * lambda1_p + lambda2_p + lambda3_p - rho_2gamma = 0.5 * rho / equations.gamma + rho_2gamma = 0.5f0 * rho / equations.gamma f1p = rho_2gamma * alpha_p f2p = rho_2gamma * alpha_p * v1 f3p = rho_2gamma * (alpha_p * v2 + a * (lambda2_p - lambda3_p)) f4p = rho_2gamma * alpha_p * v3 f5p = rho_2gamma * - (alpha_p * 0.5 * (v1^2 + v2^2 + v3^2) + a * v2 * (lambda2_p - lambda3_p) + (alpha_p * 0.5f0 * (v1^2 + v2^2 + v3^2) + + a * v2 * + (lambda2_p - lambda3_p) + a^2 * (lambda2_p + lambda3_p) * equations.inv_gamma_minus_one) else # orientation == 3 lambda1 = v3 @@ -886,13 +902,15 @@ end alpha_p = 2 * (equations.gamma - 1) * lambda1_p + lambda2_p + lambda3_p - rho_2gamma = 0.5 * rho / equations.gamma + rho_2gamma = 0.5f0 * rho / equations.gamma f1p = rho_2gamma * alpha_p f2p = rho_2gamma * alpha_p * v1 f3p = rho_2gamma * alpha_p * v2 f4p = rho_2gamma * (alpha_p * v3 + a * (lambda2_p - lambda3_p)) f5p = rho_2gamma * - (alpha_p * 0.5 * (v1^2 + v2^2 + v3^2) + a * v3 * (lambda2_p - lambda3_p) + (alpha_p * 0.5f0 * (v1^2 + v2^2 + v3^2) + + a * v3 * + (lambda2_p - lambda3_p) + a^2 * (lambda2_p + lambda3_p) * equations.inv_gamma_minus_one) end return SVector(f1p, f2p, f3p, f4p, f5p) @@ -905,7 +923,7 @@ end v2 = rho_v2 / rho v3 = rho_v3 / rho p = (equations.gamma - 1) * - (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3)) + (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3)) a = sqrt(equations.gamma * p / rho) if orientation == 1 @@ -919,13 +937,15 @@ end alpha_m = 2 * (equations.gamma - 1) * lambda1_m + lambda2_m + lambda3_m - rho_2gamma = 0.5 * rho / equations.gamma + rho_2gamma = 0.5f0 * rho / equations.gamma f1m = rho_2gamma * alpha_m f2m = rho_2gamma * (alpha_m * v1 + a * (lambda2_m - lambda3_m)) f3m = rho_2gamma * alpha_m * v2 f4m = rho_2gamma * alpha_m * v3 f5m = rho_2gamma * - (alpha_m * 0.5 * (v1^2 + v2^2 + v3^2) + a * v1 * (lambda2_m - lambda3_m) + (alpha_m * 0.5f0 * (v1^2 + v2^2 + v3^2) + + a * v1 * + (lambda2_m - lambda3_m) + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one) elseif orientation == 2 lambda1 = v2 @@ -938,13 +958,15 @@ end alpha_m = 2 * (equations.gamma - 1) * lambda1_m + lambda2_m + lambda3_m - rho_2gamma = 0.5 * rho / equations.gamma + rho_2gamma = 0.5f0 * rho / equations.gamma f1m = rho_2gamma * alpha_m f2m = rho_2gamma * alpha_m * v1 f3m = rho_2gamma * (alpha_m * v2 + a * (lambda2_m - lambda3_m)) f4m = rho_2gamma * alpha_m * v3 f5m = rho_2gamma * - (alpha_m * 0.5 * (v1^2 + v2^2 + v3^2) + a * v2 * (lambda2_m - lambda3_m) + (alpha_m * 0.5f0 * (v1^2 + v2^2 + v3^2) + + a * v2 * + (lambda2_m - lambda3_m) + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one) else # orientation == 3 lambda1 = v3 @@ -957,13 +979,15 @@ end alpha_m = 2 * (equations.gamma - 1) * lambda1_m + lambda2_m + lambda3_m - rho_2gamma = 0.5 * rho / equations.gamma + rho_2gamma = 0.5f0 * rho / equations.gamma f1m = rho_2gamma * alpha_m f2m = rho_2gamma * alpha_m * v1 f3m = rho_2gamma * alpha_m * v2 f4m = rho_2gamma * (alpha_m * v3 + a * (lambda2_m - lambda3_m)) f5m = rho_2gamma * - (alpha_m * 0.5 * (v1^2 + v2^2 + v3^2) + a * v3 * (lambda2_m - lambda3_m) + (alpha_m * 0.5f0 * (v1^2 + v2^2 + v3^2) + + a * v3 * + (lambda2_m - lambda3_m) + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one) end return SVector(f1m, f2m, f3m, f4m, f5m) @@ -1001,9 +1025,9 @@ References: v_rr = v3_rr end - rho = 0.5 * (rho_ll + rho_rr) - p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) - v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) + rho = 0.5f0 * (rho_ll + rho_rr) + p = 0.5f0 * (p_ll + p_rr) - 0.5f0 * c * rho * (v_rr - v_ll) + v = 0.5f0 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) # We treat the energy term analogous to the potential temperature term in the paper by # Chen et al., i.e. we use p_ll and p_rr, and not p @@ -1043,9 +1067,9 @@ end # and then multiplying v by `norm_` again, but this version is slightly faster. norm_ = norm(normal_direction) - rho = 0.5 * (rho_ll + rho_rr) - p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_ - v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ + rho = 0.5f0 * (rho_ll + rho_rr) + p = 0.5f0 * (p_ll + p_rr) - 0.5f0 * c * rho * (v_rr - v_ll) / norm_ + v = 0.5f0 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ # We treat the energy term analogous to the potential temperature term in the paper by # Chen et al., i.e. we use p_ll and p_rr, and not p @@ -1248,7 +1272,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, v3_ll = rho_v3_ll / rho_ll e_ll = rho_e_ll / rho_ll p_ll = (equations.gamma - 1) * - (rho_e_ll - 1 / 2 * rho_ll * (v1_ll^2 + v2_ll^2 + v3_ll^2)) + (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2 + v2_ll^2 + v3_ll^2)) c_ll = sqrt(equations.gamma * p_ll / rho_ll) v1_rr = rho_v1_rr / rho_rr @@ -1256,7 +1280,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, v3_rr = rho_v3_rr / rho_rr e_rr = rho_e_rr / rho_rr p_rr = (equations.gamma - 1) * - (rho_e_rr - 1 / 2 * rho_rr * (v1_rr^2 + v2_rr^2 + v3_rr^2)) + (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2 + v2_rr^2 + v3_rr^2)) c_rr = sqrt(equations.gamma * p_rr / rho_rr) # Obtain left and right fluxes @@ -1285,19 +1309,19 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, H_ll = (rho_e_ll + p_ll) / rho_ll H_rr = (rho_e_rr + p_rr) / rho_rr H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5f0 * vel_roe_mag)) Ssl = min(vel_L - c_ll, vel_roe - c_roe) Ssr = max(vel_R + c_rr, vel_roe + c_roe) sMu_L = Ssl - vel_L sMu_R = Ssr - vel_R - if Ssl >= 0.0 + if Ssl >= 0 f1 = f_ll[1] f2 = f_ll[2] f3 = f_ll[3] f4 = f_ll[4] f5 = f_ll[5] - elseif Ssr <= 0.0 + elseif Ssr <= 0 f1 = f_rr[1] f2 = f_rr[2] f3 = f_rr[3] @@ -1306,7 +1330,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, else SStar = (p_rr - p_ll + rho_ll * vel_L * sMu_L - rho_rr * vel_R * sMu_R) / (rho_ll * sMu_L - rho_rr * sMu_R) - if Ssl <= 0.0 <= SStar + if Ssl <= 0 <= SStar densStar = rho_ll * sMu_L / (Ssl - SStar) enerStar = e_ll + (SStar - vel_L) * (SStar + p_ll / (rho_ll * sMu_L)) UStar1 = densStar @@ -1398,20 +1422,20 @@ function flux_hllc(u_ll, u_rr, normal_direction::AbstractVector, H_rr = (u_rr[5] + p_rr) / rho_rr H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) * norm_ + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5f0 * vel_roe_mag)) * norm_ Ssl = min(v_dot_n_ll - c_ll, vel_roe - c_roe) Ssr = max(v_dot_n_rr + c_rr, vel_roe + c_roe) sMu_L = Ssl - v_dot_n_ll sMu_R = Ssr - v_dot_n_rr - if Ssl >= 0.0 + if Ssl >= 0 f1 = f_ll[1] f2 = f_ll[2] f3 = f_ll[3] f4 = f_ll[4] f5 = f_ll[5] - elseif Ssr <= 0.0 + elseif Ssr <= 0 f1 = f_rr[1] f2 = f_rr[2] f3 = f_rr[3] @@ -1420,7 +1444,7 @@ function flux_hllc(u_ll, u_rr, normal_direction::AbstractVector, else SStar = (rho_ll * v_dot_n_ll * sMu_L - rho_rr * v_dot_n_rr * sMu_R + (p_rr - p_ll) * norm_sq) / (rho_ll * sMu_L - rho_rr * sMu_R) - if Ssl <= 0.0 <= SStar + if Ssl <= 0 <= SStar densStar = rho_ll * sMu_L / (Ssl - SStar) enerStar = e_ll + (SStar - v_dot_n_ll) * @@ -1501,22 +1525,22 @@ of the numerical flux. v_roe_mag = v1_roe^2 + v2_roe^2 + v3_roe^2 H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) * inv_sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag)) + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5f0 * v_roe_mag)) # Compute convenience constant for positivity preservation, see # https://doi.org/10.1016/0021-9991(91)90211-3 - beta = sqrt(0.5 * (equations.gamma - 1) / equations.gamma) + beta = sqrt(0.5f0 * (equations.gamma - 1) / equations.gamma) # Estimate the edges of the Riemann fan (with positivity conservation) if orientation == 1 # x-direction - SsL = min(v1_roe - c_roe, v1_ll - beta * c_ll, zero(v1_roe)) - SsR = max(v1_roe + c_roe, v1_rr + beta * c_rr, zero(v1_roe)) + SsL = min(v1_roe - c_roe, v1_ll - beta * c_ll, 0) + SsR = max(v1_roe + c_roe, v1_rr + beta * c_rr, 0) elseif orientation == 2 # y-direction - SsL = min(v2_roe - c_roe, v2_ll - beta * c_ll, zero(v2_roe)) - SsR = max(v2_roe + c_roe, v2_rr + beta * c_rr, zero(v2_roe)) + SsL = min(v2_roe - c_roe, v2_ll - beta * c_ll, 0) + SsR = max(v2_roe + c_roe, v2_rr + beta * c_rr, 0) else # z-direction - SsL = min(v3_roe - c_roe, v3_ll - beta * c_ll, zero(v3_roe)) - SsR = max(v3_roe + c_roe, v3_rr + beta * c_rr, zero(v3_roe)) + SsL = min(v3_roe - c_roe, v3_ll - beta * c_ll, 0) + SsR = max(v3_roe + c_roe, v3_rr + beta * c_rr, 0) end return SsL, SsR @@ -1571,15 +1595,15 @@ of the numerical flux. v_roe_mag = v1_roe^2 + v2_roe^2 + v3_roe^2 H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) * inv_sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag)) * norm_ + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5f0 * v_roe_mag)) * norm_ # Compute convenience constant for positivity preservation, see # https://doi.org/10.1016/0021-9991(91)90211-3 - beta = sqrt(0.5 * (equations.gamma - 1) / equations.gamma) + beta = sqrt(0.5f0 * (equations.gamma - 1) / equations.gamma) # Estimate the edges of the Riemann fan (with positivity conservation) - SsL = min(v_roe - c_roe, v_dot_n_ll - beta * c_ll, zero(v_roe)) - SsR = max(v_roe + c_roe, v_dot_n_rr + beta * c_rr, zero(v_roe)) + SsL = min(v_roe - c_roe, v_dot_n_ll - beta * c_ll, 0) + SsR = max(v_roe + c_roe, v_dot_n_rr + beta * c_rr, 0) return SsL, SsR end @@ -1599,7 +1623,7 @@ end v2 = rho_v2 / rho v3 = rho_v3 / rho p = (equations.gamma - 1) * - (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3)) + (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3)) return SVector(rho, v1, v2, v3, p) end @@ -1612,11 +1636,12 @@ end v2 = rho_v2 / rho v3 = rho_v3 / rho v_square = v1^2 + v2^2 + v3^2 - p = (equations.gamma - 1) * (rho_e - 0.5 * rho * v_square) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho * v_square) s = log(p) - equations.gamma * log(rho) rho_p = rho / p - w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - 0.5 * rho_p * v_square + w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - + 0.5f0 * rho_p * v_square w2 = rho_p * v1 w3 = rho_p * v2 w4 = rho_p * v3 @@ -1658,7 +1683,7 @@ end rho_v2 = rho * v2 rho_v3 = rho * v3 rho_e = p * equations.inv_gamma_minus_one + - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e) end @@ -1669,14 +1694,14 @@ end @inline function pressure(u, equations::CompressibleEulerEquations3D) rho, rho_v1, rho_v2, rho_v3, rho_e = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho) return p end @inline function density_pressure(u, equations::CompressibleEulerEquations3D) rho, rho_v1, rho_v2, rho_v3, rho_e = u rho_times_p = (equations.gamma - 1) * - (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2)) + (rho * rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2)) return rho_times_p end @@ -1711,7 +1736,7 @@ end # Calculate kinetic energy for a conservative state `cons` @inline function energy_kinetic(u, equations::CompressibleEulerEquations3D) rho, rho_v1, rho_v2, rho_v3, _ = u - return 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + return 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho end # Calculate internal energy for a conservative state `cons` diff --git a/test/runtests.jl b/test/runtests.jl index 49f0977bb70..836488d0d8e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -108,6 +108,7 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc_part1" include("test_unit.jl") + include("test_type.jl") include("test_visualization.jl") end diff --git a/test/test_type.jl b/test/test_type.jl new file mode 100644 index 00000000000..de02ec47110 --- /dev/null +++ b/test/test_type.jl @@ -0,0 +1,502 @@ +module TestType + +using Test +using Trixi + +include("test_trixi.jl") + +# Start with a clean environment: remove Trixi.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) + +# Run unit tests for various equations +@testset "Test Type Stability" begin + @timed_testset "Acoustic Perturbation 2D" begin + for RealT in (Float32, Float64) + v_mean_global = (zero(RealT), zero(RealT)) + c_mean_global = one(RealT) + rho_mean_global = one(RealT) + equations = @inferred AcousticPerturbationEquations2D(v_mean_global, + c_mean_global, + rho_mean_global) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), + one(RealT), + one(RealT), one(RealT)) + orientations = [1, 2] + directions = [1, 2, 3, 4] + normal_direction = SVector(one(RealT), zero(RealT)) + + surface_flux_function = flux_lax_friedrichs + dissipation = DissipationLocalLaxFriedrichs() + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_gauss(x, t, equations)) == RealT + + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + + for orientation in orientations + 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 + else + @test eltype(@inferred boundary_condition_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations)) == RealT + end + end + end + + 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 + else + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + normal_direction, x, t, + surface_flux_function, + equations)) == + RealT + end + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == + RealT + @test eltype(@inferred dissipation(u_ll, u_rr, normal_direction, equations)) == + RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == + RealT + @test eltype(@inferred dissipation(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 + end + end + + @timed_testset "Compressible Euler 1D" begin + for RealT in (Float32, Float64) + # set gamma = 2 for the coupling convergence test + equations = @inferred CompressibleEulerEquations1D(RealT(2)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT)) + orientation = 1 + directions = [1, 2] + cons = SVector(one(RealT), one(RealT), one(RealT)) + + surface_flux_function = flux_lax_friedrichs + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_density_wave(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_eoc_test_coupled_euler_gravity(x, t, + equations)) == + RealT + + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + + for direction in directions + @test eltype(@inferred boundary_condition_slip_wall(u_inner, orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + end + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_shima_etal(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred flux_kennedy_gruber(u_ll, u_rr, orientation, equations)) == + RealT + @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)) == + 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)) == + RealT + else + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT + end + + @test eltype(eltype(@inferred splitting_steger_warming(u, orientation, + equations))) == + RealT + @test eltype(eltype(@inferred splitting_vanleer_haenel(u, orientation, + equations))) == + RealT + @test eltype(eltype(@inferred splitting_coirier_vanleer(u, orientation, + equations))) == + RealT + + @test eltype(@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 min_max_speed_davis(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred min_max_speed_einfeldt(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 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 + end + end + + @timed_testset "Compressible Euler 2D" begin + for RealT in (Float32, Float64) + # set gamma = 2 for the coupling convergence test\ + equations = @inferred CompressibleEulerEquations2D(RealT(2)) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT), + one(RealT)) + orientations = [1, 2] + directions = [1, 2, 3, 4] + normal_direction = SVector(one(RealT), zero(RealT)) + cons = SVector(one(RealT), one(RealT), one(RealT), one(RealT)) + + surface_flux_function = flux_lax_friedrichs + flux_lmars = FluxLMARS(340) + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_density_wave(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_eoc_test_coupled_euler_gravity(x, t, + equations)) == + RealT + + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + @test eltype(@inferred source_terms_eoc_test_coupled_euler_gravity(u, x, t, + equations)) == + RealT + @test eltype(@inferred source_terms_eoc_test_euler(u, x, t, equations)) == + RealT + + for orientation in orientations + for direction in directions + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + orientation, + direction, x, t, + surface_flux_function, + equations)) == RealT + end + end + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_shima_etal(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_kennedy_gruber(u_ll, u_rr, normal_direction, + equations)) == + RealT + @test eltype(@inferred flux_lmars(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_hllc(u_ll, u_rr, normal_direction, equations)) == + RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(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)) == + RealT + else + @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, + equations)) == RealT + end + + @test eltype(eltype(@inferred splitting_drikakis_tsangaris(u, normal_direction, + equations))) == RealT + @test eltype(eltype(@inferred splitting_vanleer_haenel(u, normal_direction, + equations))) == + RealT + @test eltype(eltype(@inferred splitting_lax_friedrichs(u, normal_direction, + equations))) == + RealT + + @test eltype(@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, + equations)) == + RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction, + equations)) == + RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, normal_direction, + equations)) == + RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_shima_etal(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred flux_kennedy_gruber(u_ll, u_rr, orientation, + equations)) == + RealT + @test eltype(@inferred flux_lmars(u_ll, u_rr, orientation, equations)) == + RealT + @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)) == + 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)) == + RealT + else + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT + end + + @test eltype(eltype(@inferred splitting_steger_warming(u, orientation, + equations))) == + RealT + @test eltype(eltype(@inferred splitting_drikakis_tsangaris(u, orientation, + equations))) == + RealT + @test eltype(eltype(@inferred splitting_vanleer_haenel(u, orientation, + equations))) == + RealT + @test eltype(eltype(@inferred splitting_lax_friedrichs(u, orientation, + equations))) == + RealT + + @test eltype(@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 min_max_speed_davis(u_ll, u_rr, orientation, + equations)) == + RealT + @test eltype(@inferred min_max_speed_einfeldt(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 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 + # TODO: test `gradient_conservative`, not necessary but good to have + end + end + + @timed_testset "Compressible Euler 3D" begin + for RealT in (Float32, Float64) + # set gamma = 2 for the coupling convergence test + equations = @inferred CompressibleEulerEquations3D(RealT(2)) + + x = SVector(zero(RealT), zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT)) + orientations = [1, 2, 3] + directions = [1, 2, 3, 4, 5, 6] + normal_direction = SVector(one(RealT), zero(RealT), zero(RealT)) + cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT)) + + surface_flux_function = flux_lax_friedrichs + flux_lmars = FluxLMARS(340) + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == 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 initial_condition_eoc_test_coupled_euler_gravity(x, t, + equations)) == + RealT + + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + @test eltype(@inferred source_terms_eoc_test_coupled_euler_gravity(u, x, t, + equations)) == + RealT + @test eltype(@inferred source_terms_eoc_test_euler(u, x, t, equations)) == RealT + + for orientation in orientations + for direction in directions + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + orientation, + direction, x, t, + surface_flux_function, + equations)) == RealT + end + end + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_shima_etal(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_kennedy_gruber(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred flux_lmars(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_hllc(u_ll, u_rr, normal_direction, equations)) == + RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(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)) == + 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, + equations)) == + RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, normal_direction, + equations)) == RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_shima_etal(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred flux_kennedy_gruber(u_ll, u_rr, orientation, + equations)) == + RealT + @test eltype(@inferred flux_lmars(u_ll, u_rr, orientation, equations)) == + RealT + @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)) == 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)) == + RealT + else + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT + end + + @test eltype(eltype(@inferred splitting_steger_warming(u, orientation, + equations))) == + RealT + + @test eltype(@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 min_max_speed_davis(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_einfeldt(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 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 + end + end +end + +end # module