From cee9c6fbeff2dd0b06fd271cd1f6882b679813e4 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Fri, 19 Apr 2024 15:11:14 -0700 Subject: [PATCH 01/21] start with src/equations --- src/equations/acoustic_perturbation_2d.jl | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/equations/acoustic_perturbation_2d.jl b/src/equations/acoustic_perturbation_2d.jl index f4ce770e1e9..f0b44b26870 100644 --- a/src/equations/acoustic_perturbation_2d.jl +++ b/src/equations/acoustic_perturbation_2d.jl @@ -116,7 +116,9 @@ function initial_condition_constant(x, t, equations::AcousticPerturbationEquatio v2_prime = 0.0 p_prime_scaled = 0.0 - return SVector(v1_prime, v2_prime, p_prime_scaled, global_mean_vars(equations)...) + return SVector(map(v -> convert(eltype(x), v), + (v1_prime, v2_prime, p_prime_scaled, + global_mean_vars(equations)...))) end """ @@ -141,7 +143,7 @@ function initial_condition_convergence_test(x, t, prim = SVector(v1_prime, v2_prime, p_prime, global_mean_vars(equations)...) - return prim2cons(prim, equations) + return map(v -> convert(eltype(x), v), prim2cons(prim, equations)) end """ @@ -170,7 +172,7 @@ function source_terms_convergence_test(u, x, t, du4 = du5 = du6 = du7 = 0.0 - return SVector(du1, du2, du3, du4, du5, du6, du7) + return SVector(map(v -> convert(eltype(u), v), (du1, du2, du3, du4, du5, du6, du7))) end """ @@ -185,7 +187,7 @@ function initial_condition_gauss(x, t, equations::AcousticPerturbationEquations2 prim = SVector(v1_prime, v2_prime, p_prime, global_mean_vars(equations)...) - return prim2cons(prim, equations) + return map(v -> convert(eltype(x), v), prim2cons(prim, equations)) end """ @@ -216,7 +218,7 @@ function boundary_condition_wall(u_inner, orientation, direction, x, t, flux = surface_flux_function(u_boundary, u_inner, orientation, equations) end - return flux + return flux # TODO: Check surface_flux_function end """ @@ -247,7 +249,7 @@ function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, # calculate the boundary flux flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations) - return flux + return map(v -> convert(eltype(u_inner), v), flux) end # Calculate 1D flux for a single point @@ -346,7 +348,8 @@ end equations) diss = -0.5 * λ * (u_rr - u_ll) z = zero(eltype(u_ll)) - return SVector(diss[1], diss[2], diss[3], z, z, z, z) + return SVector(v -> convert(eltype(u_ll), v), + (diss[1], diss[2], diss[3], z, z, z, z)) end @inline have_constant_speed(::AcousticPerturbationEquations2D) = False() From 195da8de2c2c7ad7f9872d7de8eb96f38116c1c1 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Sun, 28 Apr 2024 16:03:13 -0700 Subject: [PATCH 02/21] revise after the first review --- src/equations/acoustic_perturbation_2d.jl | 78 ++++++++++++----------- 1 file changed, 42 insertions(+), 36 deletions(-) diff --git a/src/equations/acoustic_perturbation_2d.jl b/src/equations/acoustic_perturbation_2d.jl index f0b44b26870..4e3acb4909e 100644 --- a/src/equations/acoustic_perturbation_2d.jl +++ b/src/equations/acoustic_perturbation_2d.jl @@ -112,13 +112,12 @@ 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 + RealT = eltype(x) + v1_prime = zero(RealT) + v2_prime = zero(RealT) + p_prime_scaled = zero(RealT) - return SVector(map(v -> convert(eltype(x), v), - (v1_prime, v2_prime, p_prime_scaled, - global_mean_vars(equations)...))) + return SVector(v1_prime, v2_prime, p_prime_scaled, global_mean_vars(equations)...) end """ @@ -129,12 +128,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 / 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 @@ -143,7 +143,7 @@ function initial_condition_convergence_test(x, t, prim = SVector(v1_prime, v2_prime, p_prime, global_mean_vars(equations)...) - return map(v -> convert(eltype(x), v), prim2cons(prim, equations)) + return prim2cons(prim, equations) end """ @@ -156,12 +156,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 / 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 @@ -170,9 +171,9 @@ 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 = zero(RealT) - return SVector(map(v -> convert(eltype(u), v), (du1, du2, du3, du4, du5, du6, du7))) + return SVector(du1, du2, du3, du4, du5, du6, du7) end """ @@ -181,13 +182,14 @@ 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 + RealT = eltype(x) + v1_prime = zero(RealT) + v2_prime = zero(RealT) p_prime = exp(-4 * (x[1]^2 + x[2]^2)) prim = SVector(v1_prime, v2_prime, p_prime, global_mean_vars(equations)...) - return map(v -> convert(eltype(x), v), prim2cons(prim, equations)) + return prim2cons(prim, equations) end """ @@ -218,7 +220,7 @@ function boundary_condition_wall(u_inner, orientation, direction, x, t, flux = surface_flux_function(u_boundary, u_inner, orientation, equations) end - return flux # TODO: Check surface_flux_function + return flux end """ @@ -242,14 +244,14 @@ 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 flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations) - return map(v -> convert(eltype(u_inner), v), flux) + return flux end # Calculate 1D flux for a single point @@ -259,13 +261,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 @@ -274,7 +277,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 = zero(RealT) return SVector(f1, f2, f3, f4, f5, f6, f7) end @@ -315,7 +318,8 @@ 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)) + RealT = eltype(u) + f4 = f5 = f6 = f7 = zero(RealT) return SVector(f1, f2, f3, f4, f5, f6, f7) end @@ -346,10 +350,12 @@ 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)) - return SVector(v -> convert(eltype(u_ll), v), - (diss[1], diss[2], diss[3], z, z, z, z)) + + RealT = eltype(u_ll) + diss = -convert(RealT, 0.5) * λ * (u_rr - u_ll) + z = zero(RealT) + + return SVector(diss[1], diss[2], diss[3], z, z, z, z) end @inline have_constant_speed(::AcousticPerturbationEquations2D) = False() From 2801e538c8951ae692504bdc548fe844a6872cbc Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Sun, 28 Apr 2024 16:06:28 -0700 Subject: [PATCH 03/21] format src/equations --- src/equations/acoustic_perturbation_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/acoustic_perturbation_2d.jl b/src/equations/acoustic_perturbation_2d.jl index 4e3acb4909e..fe68797f4e8 100644 --- a/src/equations/acoustic_perturbation_2d.jl +++ b/src/equations/acoustic_perturbation_2d.jl @@ -350,7 +350,7 @@ end equations::AcousticPerturbationEquations2D) λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations) - + RealT = eltype(u_ll) diss = -convert(RealT, 0.5) * λ * (u_rr - u_ll) z = zero(RealT) From e2e413a4c5801713bd46eb72ab12484d9d2611db Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Tue, 30 Apr 2024 11:17:22 -0700 Subject: [PATCH 04/21] another version --- src/equations/acoustic_perturbation_2d.jl | 28 +++++++++-------------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/src/equations/acoustic_perturbation_2d.jl b/src/equations/acoustic_perturbation_2d.jl index fe68797f4e8..07c5d176175 100644 --- a/src/equations/acoustic_perturbation_2d.jl +++ b/src/equations/acoustic_perturbation_2d.jl @@ -112,10 +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) - RealT = eltype(x) - v1_prime = zero(RealT) - v2_prime = zero(RealT) - p_prime_scaled = zero(RealT) + v1_prime = 0 + v2_prime = 0 + p_prime_scaled = 0 return SVector(v1_prime, v2_prime, p_prime_scaled, global_mean_vars(equations)...) end @@ -182,9 +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) - RealT = eltype(x) - v1_prime = zero(RealT) - v2_prime = zero(RealT) + 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)...) @@ -261,14 +259,13 @@ 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(RealT) + f2 = 0 f3 = rho_mean * v1_prime + v1_mean * p_prime_scaled else - f1 = zero(RealT) + f1 = 0 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 @@ -277,7 +274,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(RealT) + f4 = f5 = f6 = f7 = 0 return SVector(f1, f2, f3, f4, f5, f6, f7) end @@ -318,8 +315,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. - RealT = eltype(u) - f4 = f5 = f6 = f7 = zero(RealT) + f4 = f5 = f6 = f7 = 0 return SVector(f1, f2, f3, f4, f5, f6, f7) end @@ -350,10 +346,8 @@ end equations::AcousticPerturbationEquations2D) λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations) - - RealT = eltype(u_ll) - diss = -convert(RealT, 0.5) * λ * (u_rr - u_ll) - z = zero(RealT) + diss = -0.5f0 * λ * (u_rr - u_ll) + z = 0 return SVector(diss[1], diss[2], diss[3], z, z, z, z) end From ccc4a064102c5f74be830c4f33251e2b66425fb0 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 2 May 2024 10:45:49 -0700 Subject: [PATCH 05/21] revise after the second review --- src/equations/acoustic_perturbation_2d.jl | 6 +- src/equations/compressible_euler_1d.jl | 219 +++++++++++----------- 2 files changed, 115 insertions(+), 110 deletions(-) diff --git a/src/equations/acoustic_perturbation_2d.jl b/src/equations/acoustic_perturbation_2d.jl index 07c5d176175..71cc22e48e7 100644 --- a/src/equations/acoustic_perturbation_2d.jl +++ b/src/equations/acoustic_perturbation_2d.jl @@ -130,7 +130,7 @@ function initial_condition_convergence_test(x, t, RealT = eltype(x) a = 1 c = 2 - L = 2 + L = convert(RealT, 2) f = 2 / L A = convert(RealT, 0.2) omega = 2 * convert(RealT, pi) * f @@ -158,7 +158,7 @@ function source_terms_convergence_test(u, x, t, RealT = eltype(u) a = 1 c = 2 - L = 2 + L = convert(RealT, 2) f = 2 / L A = convert(RealT, 0.2) omega = 2 * convert(RealT, pi) * f @@ -170,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 = zero(RealT) + du4 = du5 = du6 = du7 = 0 return SVector(du1, du2, du3, du4, du5, du6, du7) end diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index a50c896cd1c..f4de0754ab3 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 - L = 2 + A = convert(RealT, 0.1) + L = convert(RealT, 2) f = 1 / L - ω = 2 * pi * f + ω = 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 - L = 2 + A = convert(RealT, 0.1) + L = convert(RealT, 2) f = 1 / L - ω = 2 * pi * f + ω = 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,8 +133,9 @@ 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 @@ -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.5 ? 1 : convert(RealT, 1.1691) + v1 = r > 0.5 ? 0 : convert(RealT, 0.1882) * cos_phi + p = r > 0.5 ? 1 : 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_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 + rho_avg = 0.5f0 * (rho_ll + rho_rr) + rho_mean = ln_mean(rho_ll, rho_rr) # TODO: Modify ln_mean() from math.jl + beta_mean = ln_mean(beta_ll, beta_rr) # TODO: Modify ln_mean() from math.jl + 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) @@ -388,22 +392,22 @@ See also rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations) # Compute the necessary mean values - rho_mean = ln_mean(rho_ll, rho_rr) + rho_mean = ln_mean(rho_ll, rho_rr) # TODO: Modify ln_mean() from math.jl # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)` # in exact arithmetic since # 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) + inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) # TODO: Modify inv_ln_mean() from math.jl + 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` From 388917b414d64a4e45bf92618ba64a058c2870ce Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Fri, 3 May 2024 10:48:49 -0700 Subject: [PATCH 06/21] apply suggestions from code review - comments Co-authored-by: Hendrik Ranocha --- src/equations/acoustic_perturbation_2d.jl | 4 ++-- src/equations/compressible_euler_1d.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/equations/acoustic_perturbation_2d.jl b/src/equations/acoustic_perturbation_2d.jl index 71cc22e48e7..d35ef1d2b42 100644 --- a/src/equations/acoustic_perturbation_2d.jl +++ b/src/equations/acoustic_perturbation_2d.jl @@ -130,7 +130,7 @@ function initial_condition_convergence_test(x, t, RealT = eltype(x) a = 1 c = 2 - L = convert(RealT, 2) + L = convert(RealT, 2) # since we divide by L below f = 2 / L A = convert(RealT, 0.2) omega = 2 * convert(RealT, pi) * f @@ -158,7 +158,7 @@ function source_terms_convergence_test(u, x, t, RealT = eltype(u) a = 1 c = 2 - L = convert(RealT, 2) + L = convert(RealT, 2) # since we divide by L below f = 2 / L A = convert(RealT, 0.2) omega = 2 * convert(RealT, pi) * f diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index f4de0754ab3..92e6b93ce14 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -72,7 +72,7 @@ function initial_condition_convergence_test(x, t, RealT = eltype(x) c = 2 A = convert(RealT, 0.1) - L = convert(RealT, 2) + L = convert(RealT, 2) # since we divide by L below f = 1 / L ω = 2 * convert(RealT, pi) * f ini = c + A * sin(ω * (x[1] - t)) @@ -97,7 +97,7 @@ Source terms used for convergence tests in combination with RealT = eltype(u) c = 2 A = convert(RealT, 0.1) - L = convert(RealT, 2) + L = convert(RealT, 2) # since we divide by L below f = 1 / L ω = 2 * convert(RealT, pi) * f γ = equations.gamma From 77e32d0e1052bfae617bfced6cce37b738b6a304 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Fri, 3 May 2024 12:11:46 -0700 Subject: [PATCH 07/21] remove TODO --- src/equations/compressible_euler_1d.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index 92e6b93ce14..677a661b328 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -354,8 +354,8 @@ Entropy conserving two-point flux by # Compute the necessary mean values rho_avg = 0.5f0 * (rho_ll + rho_rr) - rho_mean = ln_mean(rho_ll, rho_rr) # TODO: Modify ln_mean() from math.jl - beta_mean = ln_mean(beta_ll, beta_rr) # TODO: Modify ln_mean() from math.jl + rho_mean = ln_mean(rho_ll, rho_rr) + beta_mean = ln_mean(beta_ll, beta_rr) beta_avg = 0.5f0 * (beta_ll + beta_rr) v1_avg = 0.5f0 * (v1_ll + v1_rr) p_mean = 0.5f0 * rho_avg / beta_avg @@ -392,12 +392,12 @@ See also rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations) # Compute the necessary mean values - rho_mean = ln_mean(rho_ll, rho_rr) # TODO: Modify ln_mean() from math.jl + rho_mean = ln_mean(rho_ll, rho_rr) # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)` # in exact arithmetic since # 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) # TODO: Modify inv_ln_mean() from math.jl + inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) 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) From 8fce44d652df0dc3652646d6e8b3042b68bc9c37 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Tue, 7 May 2024 17:41:20 -0700 Subject: [PATCH 08/21] fix small errors --- src/equations/compressible_euler_1d.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index 677a661b328..4f9053778b8 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -138,7 +138,7 @@ function initial_condition_density_wave(x, t, equations::CompressibleEulerEquati 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 @@ -864,8 +864,9 @@ Compactly summarized: 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, 0) - SsR = max(v_roe + c_roe, v_rr + beta * c_rr, 0) + RealT = eltype(u_ll) + SsL = min(v_roe - c_roe, v_ll - beta * c_ll, zero(RealT)) + SsR = max(v_roe + c_roe, v_rr + beta * c_rr, zero(RealT)) return SsL, SsR end From fa90a8b7f88c8d848f299cea42bbfd6a6e8cdf8d Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Tue, 7 May 2024 19:22:25 -0700 Subject: [PATCH 09/21] complete compressible_euler --- src/equations/compressible_euler_1d.jl | 6 +- src/equations/compressible_euler_2d.jl | 445 +++++++++++++------------ src/equations/compressible_euler_3d.jl | 370 ++++++++++---------- 3 files changed, 430 insertions(+), 391 deletions(-) diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index 4f9053778b8..b30cd8794e2 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -165,9 +165,9 @@ function initial_condition_weak_blast_wave(x, t, cos_phi = x_norm > 0 ? 1 : -1 # Calculate primitive variables - rho = r > 0.5 ? 1 : convert(RealT, 1.1691) - v1 = r > 0.5 ? 0 : convert(RealT, 0.1882) * cos_phi - p = r > 0.5 ? 1 : convert(RealT, 1.245) + rho = r > 0.5f0 ? 1 : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * cos_phi + p = r > 0.5f0 ? 1 : convert(RealT, 1.245) return prim2cons(SVector(rho, v1, p), equations) end diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index fa7af285aa4..3fab27b2f39 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 - L = 2 + A = convert(RealT, 0.1) + L = convert(RealT, 2) # since we divide by L below f = 1 / L - ω = 2 * pi * f + ω = 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 - L = 2 + A = convert(RealT, 0.1) + L = convert(RealT, 2) # since we divide by L below f = 1 / L - ω = 2 * pi * f + ω = 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,10 @@ 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 + rho = r > 0.5f0 ? 1 : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * cos_phi + v2 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * sin_phi + p = r > 0.5f0 ? 1 : convert(RealT, 1.245) return prim2cons(SVector(rho, v1, v2, p), equations) end @@ -190,18 +194,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 +223,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 +254,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 +314,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 +346,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), 0) else # orientation == 2 - normal_direction = SVector(0, 1) + normal_direction = SVector(zero(RealT), 1) end # compute and return the flux using `boundary_condition_slip_wall` routine above @@ -380,7 +388,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 +442,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 +475,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 +488,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 +512,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 +543,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 +573,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 +593,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 +615,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 +668,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 +680,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 +708,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 +762,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 +776,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 +794,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 +810,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 +824,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 +842,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 +898,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 +909,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 +921,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 +935,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 +946,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 +958,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 +973,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 +985,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 +1000,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 +1012,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 +1055,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 +1093,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 +1165,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 +1195,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 +1226,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 +1249,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 +1299,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 +1326,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 +1356,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 +1377,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 +1565,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 +1596,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 +1615,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 +1689,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 +1709,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 +1783,20 @@ 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) + RealT = eltype(u_ll) 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, zero(RealT)) + SsR = max(v1_roe + c_roe, v1_rr + beta * c_rr, zero(RealT)) 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, zero(RealT)) + SsR = max(v2_roe + c_roe, v2_rr + beta * c_rr, zero(RealT)) end return SsL, SsR @@ -1836,15 +1847,16 @@ 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)) + RealT = eltype(u_ll) + SsL = min(v_roe - c_roe, v_dot_n_ll - beta * c_ll, zero(RealT)) + SsR = max(v_roe + c_roe, v_dot_n_rr + beta * c_rr, zero(RealT)) return SsL, SsR end @@ -1862,7 +1874,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 +1886,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 @@ -1915,7 +1928,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 @@ -1926,7 +1939,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 @@ -1939,12 +1952,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 @@ -1974,7 +1987,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]) @@ -2013,7 +2026,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..a41ce5516bd 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 - L = 2 + A = convert(RealT, 0.1) + L = convert(RealT, 2) # since we divide by L below f = 1 / L - ω = 2 * pi * f + ω = 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 - L = 2 + A = convert(RealT, 0.1) + L = convert(RealT, 2) # since we divide by L below f = 1 / L - ω = 2 * pi * f + ω = 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 ? 1 : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * cos(phi) * sin(theta) + v2 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * sin(phi) * sin(theta) + v3 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * cos(theta) + p = r > 0.5f0 ? 1 : 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), 0, 0) elseif orientation == 2 - normal_direction = SVector(0.0, 1.0, 0.0) + normal_direction = SVector(zero(RealT), 1, 0) else # orientation == 3 - normal_direction = SVector(0.0, 0.0, 1.0) + normal_direction = SVector(zero(RealT), 0, 1) 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 @@ -448,30 +456,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 +501,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 +516,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 +540,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 +587,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 +624,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 +646,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 +681,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 +736,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 +750,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 +758,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 +766,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 +789,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 +845,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 +859,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 +880,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 +901,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 +922,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 +936,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 +957,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 +978,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 +1024,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 +1066,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 +1271,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 +1279,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 +1308,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 +1329,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 +1421,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 +1443,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 +1524,23 @@ 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) + RealT = eltype(u_ll) 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, zero(RealT)) + SsR = max(v1_roe + c_roe, v1_rr + beta * c_rr, zero(RealT)) 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, zero(RealT)) + SsR = max(v2_roe + c_roe, v2_rr + beta * c_rr, zero(RealT)) 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, zero(RealT)) + SsR = max(v3_roe + c_roe, v3_rr + beta * c_rr, zero(RealT)) end return SsL, SsR @@ -1571,15 +1595,16 @@ 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)) + RealT = eltype(u_ll) + SsL = min(v_roe - c_roe, v_dot_n_ll - beta * c_ll, zero(RealT)) + SsR = max(v_roe + c_roe, v_dot_n_rr + beta * c_rr, zero(RealT)) return SsL, SsR end @@ -1599,7 +1624,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 +1637,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 +1684,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 +1695,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 +1737,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` From 832ad16bf9c93b84597988cd002c14bd4567445d Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Tue, 7 May 2024 19:51:27 -0700 Subject: [PATCH 10/21] fix error and format --- src/equations/compressible_euler_2d.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 3fab27b2f39..0be061d0f02 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -173,6 +173,7 @@ function initial_condition_weak_blast_wave(x, t, sin_phi, cos_phi = sincos(phi) # Calculate primitive variables + RealT = eltype(x) rho = r > 0.5f0 ? 1 : convert(RealT, 1.1691) v1 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * cos_phi v2 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * sin_phi From 08272e0465a73c3c2bf49a2e2240b415f3035194 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Wed, 8 May 2024 10:35:19 -0700 Subject: [PATCH 11/21] revise min max --- src/equations/compressible_euler_1d.jl | 5 ++--- src/equations/compressible_euler_2d.jl | 14 ++++++-------- src/equations/compressible_euler_3d.jl | 18 ++++++++---------- 3 files changed, 16 insertions(+), 21 deletions(-) diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index b30cd8794e2..9ae04f43ff3 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -864,9 +864,8 @@ Compactly summarized: beta = sqrt(0.5f0 * (equations.gamma - 1) / equations.gamma) # Estimate the edges of the Riemann fan (with positivity conservation) - RealT = eltype(u_ll) - SsL = min(v_roe - c_roe, v_ll - beta * c_ll, zero(RealT)) - SsR = max(v_roe + c_roe, v_rr + beta * c_rr, zero(RealT)) + 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 diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 0be061d0f02..44b2bac6f48 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1791,13 +1791,12 @@ of the numerical flux. beta = sqrt(0.5f0 * (equations.gamma - 1) / equations.gamma) # Estimate the edges of the Riemann fan (with positivity conservation) - RealT = eltype(u_ll) if orientation == 1 # x-direction - SsL = min(v1_roe - c_roe, v1_ll - beta * c_ll, zero(RealT)) - SsR = max(v1_roe + c_roe, v1_rr + beta * c_rr, zero(RealT)) + 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(RealT)) - SsR = max(v2_roe + c_roe, v2_rr + beta * c_rr, zero(RealT)) + 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 @@ -1855,9 +1854,8 @@ of the numerical flux. beta = sqrt(0.5f0 * (equations.gamma - 1) / equations.gamma) # Estimate the edges of the Riemann fan (with positivity conservation) - RealT = eltype(u_ll) - SsL = min(v_roe - c_roe, v_dot_n_ll - beta * c_ll, zero(RealT)) - SsR = max(v_roe + c_roe, v_dot_n_rr + beta * c_rr, zero(RealT)) + 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 diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index a41ce5516bd..e26280e091d 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -1531,16 +1531,15 @@ of the numerical flux. beta = sqrt(0.5f0 * (equations.gamma - 1) / equations.gamma) # Estimate the edges of the Riemann fan (with positivity conservation) - RealT = eltype(u_ll) if orientation == 1 # x-direction - SsL = min(v1_roe - c_roe, v1_ll - beta * c_ll, zero(RealT)) - SsR = max(v1_roe + c_roe, v1_rr + beta * c_rr, zero(RealT)) + 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(RealT)) - SsR = max(v2_roe + c_roe, v2_rr + beta * c_rr, zero(RealT)) + 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(RealT)) - SsR = max(v3_roe + c_roe, v3_rr + beta * c_rr, zero(RealT)) + 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 @@ -1602,9 +1601,8 @@ of the numerical flux. beta = sqrt(0.5f0 * (equations.gamma - 1) / equations.gamma) # Estimate the edges of the Riemann fan (with positivity conservation) - RealT = eltype(u_ll) - SsL = min(v_roe - c_roe, v_dot_n_ll - beta * c_ll, zero(RealT)) - SsR = max(v_roe + c_roe, v_dot_n_rr + beta * c_rr, zero(RealT)) + 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 From 7b548904a0677bd07cd938d6fb6d12dec6ceb6d0 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Tue, 14 May 2024 14:12:22 -0700 Subject: [PATCH 12/21] revise based on new docs --- src/equations/compressible_euler_1d.jl | 6 +++--- src/equations/compressible_euler_2d.jl | 12 ++++++------ src/equations/compressible_euler_3d.jl | 16 ++++++++-------- 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index 9ae04f43ff3..20f3ddf5f16 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -165,9 +165,9 @@ function initial_condition_weak_blast_wave(x, t, cos_phi = x_norm > 0 ? 1 : -1 # Calculate primitive variables - rho = r > 0.5f0 ? 1 : convert(RealT, 1.1691) - v1 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * cos_phi - p = r > 0.5f0 ? 1 : convert(RealT, 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 diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 4e4fc6d4b1f..53652d9d59d 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -174,10 +174,10 @@ function initial_condition_weak_blast_wave(x, t, # Calculate primitive variables RealT = eltype(x) - rho = r > 0.5f0 ? 1 : convert(RealT, 1.1691) - v1 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * cos_phi - v2 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * sin_phi - p = r > 0.5f0 ? 1 : convert(RealT, 1.245) + 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 @@ -349,9 +349,9 @@ Should be used together with [`TreeMesh`](@ref). # get the appropriate normal vector from the orientation RealT = eltype(u_inner) if orientation == 1 - normal_direction = SVector(one(RealT), 0) + normal_direction = SVector(one(RealT), zero(RealT)) else # orientation == 2 - normal_direction = SVector(zero(RealT), 1) + normal_direction = SVector(zero(RealT), one(RealT)) end # compute and return the flux using `boundary_condition_slip_wall` routine above diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index e26280e091d..8732f0fe7a4 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -157,11 +157,11 @@ function initial_condition_weak_blast_wave(x, t, theta = iszero(r) ? zero(RealT) : acos(z_norm / r) # Calculate primitive variables - rho = r > 0.5f0 ? 1 : convert(RealT, 1.1691) - v1 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * cos(phi) * sin(theta) - v2 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * sin(phi) * sin(theta) - v3 = r > 0.5f0 ? 0 : convert(RealT, 0.1882) * cos(theta) - p = r > 0.5f0 ? 1 : convert(RealT, 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 @@ -351,11 +351,11 @@ Should be used together with [`TreeMesh`](@ref). # get the appropriate normal vector from the orientation RealT = eltype(u_inner) if orientation == 1 - normal_direction = SVector(one(RealT), 0, 0) + normal_direction = SVector(one(RealT), zero(RealT), zero(RealT)) elseif orientation == 2 - normal_direction = SVector(zero(RealT), 1, 0) + normal_direction = SVector(zero(RealT), one(RealT), zero(RealT)) else # orientation == 3 - normal_direction = SVector(zero(RealT), 0, 1) + normal_direction = SVector(zero(RealT), zero(RealT), one(RealT)) end # compute and return the flux using `boundary_condition_slip_wall` routine above From c0ce534272ed686137f19cdc3bb226d2bb19c84d Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Wed, 15 May 2024 10:01:34 -0700 Subject: [PATCH 13/21] revise based on new comments --- src/equations/acoustic_perturbation_2d.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/equations/acoustic_perturbation_2d.jl b/src/equations/acoustic_perturbation_2d.jl index d35ef1d2b42..0ae3e7ecb01 100644 --- a/src/equations/acoustic_perturbation_2d.jl +++ b/src/equations/acoustic_perturbation_2d.jl @@ -259,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 = 0 + f2 = zero(RealT) f3 = rho_mean * v1_prime + v1_mean * p_prime_scaled else - f1 = 0 + 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 From 523dac0c3f4278829c428c8a1bf4a787fd3fe7d9 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 20 May 2024 18:06:09 -0700 Subject: [PATCH 14/21] start sample test --- test/test_type.jl | 48 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 test/test_type.jl diff --git a/test/test_type.jl b/test/test_type.jl new file mode 100644 index 00000000000..9deb950f42f --- /dev/null +++ b/test/test_type.jl @@ -0,0 +1,48 @@ +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 "Float32 Type Stability" begin + @timed_testset "Type stability for acoustic perturbation 2D" begin + v_mean_global = (0.0f0, 0.0f0) + c_mean_global = 1.0f0 + rho_mean_global = 1.0f0 + equations = AcousticPerturbationEquations2D(v_mean_global, c_mean_global, + rho_mean_global) + + x = SVector(0.0f0, 0.0f0) + t = 0.0f0 + u = u_ll = u_rr = SVector(1.0f0, 1.0f0, 1.0f0, 1.0f0, 1.0f0, 1.0f0, 1.0f0) + orientations = [1, 2] + normal_direction = SVector(1.0f0, 1.0f0) + + @test eltype(initial_condition_constant(x, t, equations)) == Float32 + @test eltype(initial_condition_convergence_test(x, t, equations)) == Float32 + @test eltype(initial_condition_gauss(x, t, equations)) == Float32 + + @test eltype(source_terms_convergence_test(u, x, t, equations)) == Float32 + @test eltype(source_terms_convergence_test(u, x, t, equations)) == Float32 + + for orientation in orientations + @test eltype(flux(u, orientation, equations)) == Float32 + end + + @test eltype(flux(u, normal_direction, equations)) == Float32 + + dissipation = DissipationLocalLaxFriedrichs() + # test orientation + @test eltype(dissipation(u_ll, u_rr, orientations[1], equations)) == Float32 + # test normal direction + @test eltype(dissipation(u_ll, u_rr, normal_direction, equations)) == Float32 + end +end + +end # module From 4c8fa0e4b91d835078d05890f7f3c090e679c701 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 20 May 2024 21:15:44 -0700 Subject: [PATCH 15/21] add more tests --- test/test_type.jl | 52 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 51 insertions(+), 1 deletion(-) diff --git a/test/test_type.jl b/test/test_type.jl index 9deb950f42f..044592c7556 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -11,7 +11,7 @@ isdir(outdir) && rm(outdir, recursive = true) # Run unit tests for various equations @testset "Float32 Type Stability" begin - @timed_testset "Type stability for acoustic perturbation 2D" begin + @timed_testset "Acoustic Perturbation 2D" begin v_mean_global = (0.0f0, 0.0f0) c_mean_global = 1.0f0 rho_mean_global = 1.0f0 @@ -43,6 +43,56 @@ isdir(outdir) && rm(outdir, recursive = true) # test normal direction @test eltype(dissipation(u_ll, u_rr, normal_direction, equations)) == Float32 end + + @timed_testset "Compressible Euler 1D" begin + # set to 2.0f0 for coupling convergence test + equations = CompressibleEulerEquations1D(2.0f0) + + x = SVector(0.0f0) + t = 0.0f0 + u = u_ll = u_rr = u_inner = SVector(1.0f0, 1.0f0, 1.0f0) + orientation = 1 + directions = [1, 2] # even and odd + + @test eltype(initial_condition_constant(x, t, equations)) == Float32 + @test eltype(initial_condition_convergence_test(x, t, equations)) == Float32 + @test eltype(initial_condition_density_wave(x, t, equations)) == Float32 + @test eltype(initial_condition_weak_blast_wave(x, t, equations)) == Float32 + @test eltype(initial_condition_eoc_test_coupled_euler_gravity(x, t, equations)) == + Float32 + + @test eltype(source_terms_convergence_test(u, x, t, equations)) == Float32 + + # set a surface flux function as ? + #= for direction in directions + @test eltype(boundary_condition_slip_wall(u_inner, orientation, direction, x, t, + surface_flux_function, equations)) == Float32 + end =# + + @test eltype(flux(u, orientation, equations)) == Float32 + @test eltype(flux_shima_etal(u_ll, u_rr, orientation, equations)) == Float32 + @test eltype(flux_kennedy_gruber(u_ll, u_rr, orientation, equations)) == Float32 + @test_broken eltype(flux_chandrashekar(u_ll, u_rr, orientation, equations)) == + Float32 # ln_mean + @test_broken eltype(flux_ranocha(u_ll, u_rr, orientation, equations)) == Float32 # ln_mean, inv_ln_mean + @test eltype(flux_hllc(u_ll, u_rr, orientation, equations)) == Float32 + + @test eltype(eltype(splitting_steger_warming(u, orientation, equations))) == Float32 + @test eltype(eltype(splitting_vanleer_haenel(u, orientation, equations))) == Float32 + @test eltype(eltype(splitting_coirier_vanleer(u, orientation, equations))) == + Float32 + + @test eltype(max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == Float32 + @test eltype(min_max_speed_einfeldt(u_ll, u_rr, orientation, equations)) == Float32 + @test eltype(max_abs_speeds(u, equations)) == Float32 + + @test eltype(cons2prim(u, equations)) == Float32 + @test eltype(prim2cons(u, equations)) == Float32 + @test eltype(cons2entropy(u, equations)) == Float32 + @test eltype(entropy2cons(u, equations)) == Float32 + + # more minor funciton tests... (overwork today and stop here) + end end end # module From b35595133c622db8169234e7ba59ee88e6eec430 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 20 May 2024 21:19:48 -0700 Subject: [PATCH 16/21] fix typos and comments --- test/test_type.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_type.jl b/test/test_type.jl index 044592c7556..9ef1ebe4b1e 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -84,14 +84,14 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == Float32 @test eltype(min_max_speed_einfeldt(u_ll, u_rr, orientation, equations)) == Float32 - @test eltype(max_abs_speeds(u, equations)) == Float32 + # @test eltype(max_abs_speeds(u, equations)) == Float32 # not defined ? @test eltype(cons2prim(u, equations)) == Float32 @test eltype(prim2cons(u, equations)) == Float32 @test eltype(cons2entropy(u, equations)) == Float32 @test eltype(entropy2cons(u, equations)) == Float32 - # more minor funciton tests... (overwork today and stop here) + # more minor function tests... (overwork today and stop here) end end From a5b462df6c4d60831094a71acc894e0d4e08e239 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Sat, 25 May 2024 11:13:17 -0700 Subject: [PATCH 17/21] change code review --- src/equations/acoustic_perturbation_2d.jl | 8 ++++---- src/equations/compressible_euler_1d.jl | 8 ++++---- src/equations/compressible_euler_2d.jl | 8 ++++---- src/equations/compressible_euler_3d.jl | 8 ++++---- 4 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/equations/acoustic_perturbation_2d.jl b/src/equations/acoustic_perturbation_2d.jl index 0ae3e7ecb01..1bde94a648a 100644 --- a/src/equations/acoustic_perturbation_2d.jl +++ b/src/equations/acoustic_perturbation_2d.jl @@ -130,8 +130,8 @@ function initial_condition_convergence_test(x, t, RealT = eltype(x) a = 1 c = 2 - L = convert(RealT, 2) # since we divide by L below - f = 2 / L + 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)) @@ -158,8 +158,8 @@ function source_terms_convergence_test(u, x, t, RealT = eltype(u) a = 1 c = 2 - L = convert(RealT, 2) # since we divide by L below - f = 2 / L + L = 2 + f = 2.0f0 / L A = convert(RealT, 0.2) omega = 2 * convert(RealT, pi) * f diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index 20f3ddf5f16..a71750ff98c 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -72,8 +72,8 @@ function initial_condition_convergence_test(x, t, RealT = eltype(x) c = 2 A = convert(RealT, 0.1) - L = convert(RealT, 2) # since we divide by L below - f = 1 / L + L = 2 + f = 1.0f0 / L ω = 2 * convert(RealT, pi) * f ini = c + A * sin(ω * (x[1] - t)) @@ -97,8 +97,8 @@ Source terms used for convergence tests in combination with RealT = eltype(u) c = 2 A = convert(RealT, 0.1) - L = convert(RealT, 2) # since we divide by L below - f = 1 / L + L = 2 + f = 1.0f0 / L ω = 2 * convert(RealT, pi) * f γ = equations.gamma diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 394a5d32b74..6a43416b777 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -80,8 +80,8 @@ function initial_condition_convergence_test(x, t, RealT = eltype(x) c = 2 A = convert(RealT, 0.1) - L = convert(RealT, 2) # since we divide by L below - f = 1 / L + L = 2 + f = 1.0f0 / L ω = 2 * convert(RealT, pi) * f ini = c + A * sin(ω * (x[1] + x[2] - t)) @@ -106,8 +106,8 @@ Source terms used for convergence tests in combination with RealT = eltype(u) c = 2 A = convert(RealT, 0.1) - L = convert(RealT, 2) # since we divide by L below - f = 1 / L + L = 2 + f = 1.0f0 / L ω = 2 * convert(RealT, pi) * f γ = equations.gamma diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index 8732f0fe7a4..c8476728f07 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -87,8 +87,8 @@ function initial_condition_convergence_test(x, t, RealT = eltype(x) c = 2 A = convert(RealT, 0.1) - L = convert(RealT, 2) # since we divide by L below - f = 1 / L + L = 2 + f = 1.0f0 / L ω = 2 * convert(RealT, pi) * f ini = c + A * sin(ω * (x[1] + x[2] + x[3] - t)) @@ -113,8 +113,8 @@ Source terms used for convergence tests in combination with RealT = eltype(u) c = 2 A = convert(RealT, 0.1) - L = convert(RealT, 2) # since we divide by L below - f = 1 / L + L = 2 + f = 1.0f0 / L ω = 2 * convert(RealT, pi) * f γ = equations.gamma From f43d8d42eec57301573386acd06ab76a98043087 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Sat, 25 May 2024 16:01:30 -0700 Subject: [PATCH 18/21] add to CI --- src/Trixi.jl | 2 +- src/equations/compressible_euler_2d.jl | 6 +- test/runtests.jl | 1 + test/test_type.jl | 353 +++++++++++++++++++------ 4 files changed, 282 insertions(+), 80 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index f3977f1f058..ec4ff8a3a59 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -181,7 +181,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal, hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, - FluxLaxFriedrichs, max_abs_speed_naive, + FluxLaxFriedrichs, max_abs_speed_naive, max_abs_speeds, FluxHLL, min_max_speed_naive, min_max_speed_davis, min_max_speed_einfeldt, FluxLMARS, FluxRotated, diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 6a43416b777..e964d9c2b21 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1907,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 @@ -2046,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 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 index 9ef1ebe4b1e..cad74397576 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -10,88 +10,289 @@ outdir = "out" isdir(outdir) && rm(outdir, recursive = true) # Run unit tests for various equations -@testset "Float32 Type Stability" begin +@testset "Test Type Stability" begin @timed_testset "Acoustic Perturbation 2D" begin - v_mean_global = (0.0f0, 0.0f0) - c_mean_global = 1.0f0 - rho_mean_global = 1.0f0 - equations = AcousticPerturbationEquations2D(v_mean_global, c_mean_global, - rho_mean_global) - - x = SVector(0.0f0, 0.0f0) - t = 0.0f0 - u = u_ll = u_rr = SVector(1.0f0, 1.0f0, 1.0f0, 1.0f0, 1.0f0, 1.0f0, 1.0f0) - orientations = [1, 2] - normal_direction = SVector(1.0f0, 1.0f0) - - @test eltype(initial_condition_constant(x, t, equations)) == Float32 - @test eltype(initial_condition_convergence_test(x, t, equations)) == Float32 - @test eltype(initial_condition_gauss(x, t, equations)) == Float32 - - @test eltype(source_terms_convergence_test(u, x, t, equations)) == Float32 - @test eltype(source_terms_convergence_test(u, x, t, equations)) == Float32 - - for orientation in orientations - @test eltype(flux(u, orientation, equations)) == Float32 - end + for RealT in (Float32, Float64) + v_mean_global = (zero(RealT), zero(RealT)) + c_mean_global = one(RealT) + rho_mean_global = one(RealT) + equations = 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(initial_condition_constant(x, t, equations)) == RealT + @test eltype(initial_condition_convergence_test(x, t, equations)) == RealT + @test eltype(initial_condition_gauss(x, t, equations)) == RealT + + @test eltype(source_terms_convergence_test(u, x, t, equations)) == RealT + @test eltype(source_terms_convergence_test(u, x, t, equations)) == RealT - @test eltype(flux(u, normal_direction, equations)) == Float32 + 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(boundary_condition_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations)) == RealT + end + end + end - dissipation = DissipationLocalLaxFriedrichs() - # test orientation - @test eltype(dissipation(u_ll, u_rr, orientations[1], equations)) == Float32 - # test normal direction - @test eltype(dissipation(u_ll, u_rr, normal_direction, equations)) == Float32 + 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(boundary_condition_slip_wall(u_inner, normal_direction, x, t, + surface_flux_function, equations)) == + RealT + end + + @test eltype(flux(u, normal_direction, equations)) == RealT + @test eltype(max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(dissipation(u_ll, u_rr, normal_direction, equations)) == RealT + + for orientation in orientations + @test eltype(flux(u, orientation, equations)) == RealT + @test eltype(max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(dissipation(u_ll, u_rr, orientation, equations)) == RealT + end + + @test eltype(max_abs_speeds(u, equations)) == RealT + @test eltype(cons2prim(u, equations)) == RealT + @test eltype(prim2cons(u, equations)) == RealT + @test eltype(cons2entropy(u, equations)) == RealT + end end @timed_testset "Compressible Euler 1D" begin - # set to 2.0f0 for coupling convergence test - equations = CompressibleEulerEquations1D(2.0f0) - - x = SVector(0.0f0) - t = 0.0f0 - u = u_ll = u_rr = u_inner = SVector(1.0f0, 1.0f0, 1.0f0) - orientation = 1 - directions = [1, 2] # even and odd - - @test eltype(initial_condition_constant(x, t, equations)) == Float32 - @test eltype(initial_condition_convergence_test(x, t, equations)) == Float32 - @test eltype(initial_condition_density_wave(x, t, equations)) == Float32 - @test eltype(initial_condition_weak_blast_wave(x, t, equations)) == Float32 - @test eltype(initial_condition_eoc_test_coupled_euler_gravity(x, t, equations)) == - Float32 - - @test eltype(source_terms_convergence_test(u, x, t, equations)) == Float32 - - # set a surface flux function as ? - #= for direction in directions - @test eltype(boundary_condition_slip_wall(u_inner, orientation, direction, x, t, - surface_flux_function, equations)) == Float32 - end =# - - @test eltype(flux(u, orientation, equations)) == Float32 - @test eltype(flux_shima_etal(u_ll, u_rr, orientation, equations)) == Float32 - @test eltype(flux_kennedy_gruber(u_ll, u_rr, orientation, equations)) == Float32 - @test_broken eltype(flux_chandrashekar(u_ll, u_rr, orientation, equations)) == - Float32 # ln_mean - @test_broken eltype(flux_ranocha(u_ll, u_rr, orientation, equations)) == Float32 # ln_mean, inv_ln_mean - @test eltype(flux_hllc(u_ll, u_rr, orientation, equations)) == Float32 - - @test eltype(eltype(splitting_steger_warming(u, orientation, equations))) == Float32 - @test eltype(eltype(splitting_vanleer_haenel(u, orientation, equations))) == Float32 - @test eltype(eltype(splitting_coirier_vanleer(u, orientation, equations))) == - Float32 - - @test eltype(max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == Float32 - @test eltype(min_max_speed_einfeldt(u_ll, u_rr, orientation, equations)) == Float32 - # @test eltype(max_abs_speeds(u, equations)) == Float32 # not defined ? - - @test eltype(cons2prim(u, equations)) == Float32 - @test eltype(prim2cons(u, equations)) == Float32 - @test eltype(cons2entropy(u, equations)) == Float32 - @test eltype(entropy2cons(u, equations)) == Float32 - - # more minor function tests... (overwork today and stop here) + for RealT in (Float32, Float64) + # set gamma = 2 for the coupling convergence test + equations = 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(initial_condition_constant(x, t, equations)) == RealT + @test eltype(initial_condition_convergence_test(x, t, equations)) == RealT + @test eltype(initial_condition_density_wave(x, t, equations)) == RealT + @test eltype(initial_condition_weak_blast_wave(x, t, equations)) == RealT + @test eltype(initial_condition_eoc_test_coupled_euler_gravity(x, t, equations)) == + RealT + + @test eltype(source_terms_convergence_test(u, x, t, equations)) == RealT + + for direction in directions + @test eltype(boundary_condition_slip_wall(u_inner, orientation, direction, + x, t, + surface_flux_function, equations)) == + RealT + end + + @test eltype(flux(u, orientation, equations)) == RealT + @test eltype(flux_shima_etal(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(flux_kennedy_gruber(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(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(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(flux_ranocha(u_ll, u_rr, orientation, equations)) == RealT + end + + @test eltype(eltype(splitting_steger_warming(u, orientation, equations))) == + RealT + @test eltype(eltype(splitting_vanleer_haenel(u, orientation, equations))) == + RealT + @test eltype(eltype(splitting_coirier_vanleer(u, orientation, equations))) == + RealT + + @test eltype(max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(min_max_speed_naive(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(min_max_speed_davis(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(min_max_speed_einfeldt(u_ll, u_rr, orientation, equations)) == + RealT + + @test eltype(max_abs_speeds(u, equations)) == RealT + @test eltype(cons2prim(u, equations)) == RealT + @test eltype(prim2cons(u, equations)) == RealT + @test eltype(cons2entropy(u, equations)) == RealT + @test eltype(entropy2cons(u, equations)) == RealT + @test eltype(density(u, equations)) == RealT + @test eltype(pressure(u, equations)) == RealT + @test eltype(density_pressure(u, equations)) == RealT + @test eltype(entropy(cons, equations)) == RealT + @test eltype(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 = 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), 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 + flux_lmars = FluxLMARS(340) + + @test eltype(initial_condition_constant(x, t, equations)) == RealT + @test eltype(initial_condition_convergence_test(x, t, equations)) == RealT + @test eltype(initial_condition_density_wave(x, t, equations)) == RealT + @test eltype(initial_condition_weak_blast_wave(x, t, equations)) == RealT + @test eltype(initial_condition_eoc_test_coupled_euler_gravity(x, t, equations)) == + RealT + + @test eltype(source_terms_convergence_test(u, x, t, equations)) == RealT + @test eltype(source_terms_eoc_test_coupled_euler_gravity(u, x, t, equations)) == + RealT + @test eltype(source_terms_eoc_test_euler(u, x, t, equations)) == + RealT + + for orientation in orientations + for direction in directions + @test eltype(boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations)) == RealT + end + end + + @test eltype(flux(u, normal_direction, equations)) == RealT + @test eltype(flux_shima_etal(u_ll, u_rr, normal_direction, equations)) == RealT + @test eltype(flux_kennedy_gruber(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(flux_lmars(u_ll, u_rr, normal_direction, equations)) == RealT + @test eltype(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(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(flux_ranocha(u_ll, u_rr, normal_direction, equations)) == RealT + end + + @test eltype(eltype(splitting_drikakis_tsangaris(u, normal_direction, + equations))) == RealT + @test eltype(eltype(splitting_vanleer_haenel(u, normal_direction, equations))) == + RealT + @test eltype(eltype(splitting_lax_friedrichs(u, normal_direction, equations))) == + RealT + + @test eltype(max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(min_max_speed_naive(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(min_max_speed_davis(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(min_max_speed_einfeldt(u_ll, u_rr, normal_direction, equations)) == + RealT + + for orientation in orientations + @test eltype(flux(u, orientation, equations)) == RealT + @test eltype(flux_shima_etal(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(flux_kennedy_gruber(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(flux_lmars(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(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(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(flux_ranocha(u_ll, u_rr, orientation, equations)) == RealT + end + + @test eltype(eltype(splitting_steger_warming(u, orientation, equations))) == + RealT + @test eltype(eltype(splitting_drikakis_tsangaris(u, orientation, equations))) == + RealT + @test eltype(eltype(splitting_vanleer_haenel(u, orientation, equations))) == + RealT + @test eltype(eltype(splitting_lax_friedrichs(u, orientation, equations))) == + RealT + + @test eltype(max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(min_max_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(min_max_speed_davis(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(min_max_speed_einfeldt(u_ll, u_rr, orientation, equations)) == + RealT + end + + @test eltype(max_abs_speeds(u, equations)) == RealT + @test eltype(cons2prim(u, equations)) == RealT + @test eltype(prim2cons(u, equations)) == RealT + @test eltype(cons2entropy(u, equations)) == RealT + @test eltype(entropy2cons(u, equations)) == RealT + # @test eltype(cons2entropy_guermond_etal(u, equations)) == RealT + @test eltype(density(u, equations)) == RealT + @test eltype(pressure(u, equations)) == RealT + @test eltype(density_pressure(u, equations)) == RealT + end end end From 4c499fc8efd9cce243cbe1da8971b68145499230 Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Sun, 26 May 2024 21:42:06 -0700 Subject: [PATCH 19/21] apply inferred macro Co-authored-by: Hendrik Ranocha --- test/test_type.jl | 197 +++++++++++++++++++++++----------------------- 1 file changed, 99 insertions(+), 98 deletions(-) diff --git a/test/test_type.jl b/test/test_type.jl index cad74397576..7e3c0519f89 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -16,8 +16,9 @@ isdir(outdir) && rm(outdir, recursive = true) v_mean_global = (zero(RealT), zero(RealT)) c_mean_global = one(RealT) rho_mean_global = one(RealT) - equations = AcousticPerturbationEquations2D(v_mean_global, c_mean_global, - rho_mean_global) + equations = @inferred AcousticPerturbationEquations2D(v_mean_global, + c_mean_global, + rho_mean_global) x = SVector(zero(RealT), zero(RealT)) t = zero(RealT) @@ -32,12 +33,12 @@ isdir(outdir) && rm(outdir, recursive = true) surface_flux_function = flux_lax_friedrichs dissipation = DissipationLocalLaxFriedrichs() - @test eltype(initial_condition_constant(x, t, equations)) == RealT - @test eltype(initial_condition_convergence_test(x, t, equations)) == RealT - @test eltype(initial_condition_gauss(x, t, equations)) == RealT + @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(source_terms_convergence_test(u, x, t, equations)) == RealT - @test eltype(source_terms_convergence_test(u, 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 @@ -48,7 +49,7 @@ isdir(outdir) && rm(outdir, recursive = true) surface_flux_function, equations)) == RealT else - @test eltype(boundary_condition_wall(u_inner, orientation, + @test eltype(@inferred boundary_condition_wall(u_inner, orientation, direction, x, t, surface_flux_function, equations)) == RealT @@ -63,34 +64,34 @@ isdir(outdir) && rm(outdir, recursive = true) surface_flux_function, equations)) == RealT else - @test eltype(boundary_condition_slip_wall(u_inner, normal_direction, x, t, + @test eltype(@inferred boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, equations)) == RealT end - @test eltype(flux(u, normal_direction, equations)) == RealT - @test eltype(max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == + @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(dissipation(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(flux(u, orientation, equations)) == RealT - @test eltype(max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT - @test eltype(dissipation(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(@inferred dissipation(u_ll, u_rr, orientation, equations)) == RealT end - @test eltype(max_abs_speeds(u, equations)) == RealT - @test eltype(cons2prim(u, equations)) == RealT - @test eltype(prim2cons(u, equations)) == RealT - @test eltype(cons2entropy(u, equations)) == RealT + @test eltype(@inferred 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 = CompressibleEulerEquations1D(RealT(2)) + equations = @inferred CompressibleEulerEquations1D(RealT(2)) x = SVector(zero(RealT)) t = zero(RealT) @@ -101,32 +102,32 @@ isdir(outdir) && rm(outdir, recursive = true) surface_flux_function = flux_lax_friedrichs - @test eltype(initial_condition_constant(x, t, equations)) == RealT - @test eltype(initial_condition_convergence_test(x, t, equations)) == RealT - @test eltype(initial_condition_density_wave(x, t, equations)) == RealT - @test eltype(initial_condition_weak_blast_wave(x, t, equations)) == RealT - @test eltype(initial_condition_eoc_test_coupled_euler_gravity(x, t, equations)) == + @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(source_terms_convergence_test(u, x, t, equations)) == RealT + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == RealT for direction in directions - @test eltype(boundary_condition_slip_wall(u_inner, orientation, direction, + @test eltype(@inferred boundary_condition_slip_wall(u_inner, orientation, direction, x, t, surface_flux_function, equations)) == RealT end - @test eltype(flux(u, orientation, equations)) == RealT - @test eltype(flux_shima_etal(u_ll, u_rr, orientation, equations)) == RealT - @test eltype(flux_kennedy_gruber(u_ll, u_rr, orientation, equations)) == RealT - @test eltype(flux_hllc(u_ll, u_rr, orientation, equations)) == RealT + @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(flux_chandrashekar(u_ll, u_rr, orientation, equations)) == + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, equations)) == RealT end if RealT == Float32 @@ -134,39 +135,39 @@ isdir(outdir) && rm(outdir, recursive = true) @test_broken eltype(flux_ranocha(u_ll, u_rr, orientation, equations)) == RealT else - @test eltype(flux_ranocha(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == RealT end - @test eltype(eltype(splitting_steger_warming(u, orientation, equations))) == + @test eltype(eltype(@inferred splitting_steger_warming(u, orientation, equations))) == RealT - @test eltype(eltype(splitting_vanleer_haenel(u, orientation, equations))) == + @test eltype(eltype(@inferred splitting_vanleer_haenel(u, orientation, equations))) == RealT - @test eltype(eltype(splitting_coirier_vanleer(u, orientation, equations))) == + @test eltype(eltype(@inferred splitting_coirier_vanleer(u, orientation, equations))) == RealT - @test eltype(max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT - @test eltype(min_max_speed_naive(u_ll, u_rr, orientation, equations)) == RealT - @test eltype(min_max_speed_davis(u_ll, u_rr, orientation, equations)) == RealT - @test eltype(min_max_speed_einfeldt(u_ll, u_rr, orientation, equations)) == + @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(max_abs_speeds(u, equations)) == RealT - @test eltype(cons2prim(u, equations)) == RealT - @test eltype(prim2cons(u, equations)) == RealT - @test eltype(cons2entropy(u, equations)) == RealT - @test eltype(entropy2cons(u, equations)) == RealT - @test eltype(density(u, equations)) == RealT - @test eltype(pressure(u, equations)) == RealT - @test eltype(density_pressure(u, equations)) == RealT - @test eltype(entropy(cons, equations)) == RealT - @test eltype(energy_internal(cons, equations)) == RealT + @test eltype(@inferred 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 = CompressibleEulerEquations2D(RealT(2)) + equations = @inferred CompressibleEulerEquations2D(RealT(2)) x = SVector(zero(RealT), zero(RealT)) t = zero(RealT) @@ -179,41 +180,41 @@ isdir(outdir) && rm(outdir, recursive = true) surface_flux_function = flux_lax_friedrichs flux_lmars = FluxLMARS(340) - @test eltype(initial_condition_constant(x, t, equations)) == RealT - @test eltype(initial_condition_convergence_test(x, t, equations)) == RealT - @test eltype(initial_condition_density_wave(x, t, equations)) == RealT - @test eltype(initial_condition_weak_blast_wave(x, t, equations)) == RealT - @test eltype(initial_condition_eoc_test_coupled_euler_gravity(x, t, equations)) == + @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(source_terms_convergence_test(u, x, t, equations)) == RealT - @test eltype(source_terms_eoc_test_coupled_euler_gravity(u, x, t, equations)) == + @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(source_terms_eoc_test_euler(u, x, t, equations)) == + @test eltype(@inferred source_terms_eoc_test_euler(u, x, t, equations)) == RealT for orientation in orientations for direction in directions - @test eltype(boundary_condition_slip_wall(u_inner, orientation, + @test eltype(@inferred boundary_condition_slip_wall(u_inner, orientation, direction, x, t, surface_flux_function, equations)) == RealT end end - @test eltype(flux(u, normal_direction, equations)) == RealT - @test eltype(flux_shima_etal(u_ll, u_rr, normal_direction, equations)) == RealT - @test eltype(flux_kennedy_gruber(u_ll, u_rr, normal_direction, equations)) == + @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(flux_lmars(u_ll, u_rr, normal_direction, equations)) == RealT - @test eltype(flux_hllc(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(flux_chandrashekar(u_ll, u_rr, normal_direction, equations)) == + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, normal_direction, equations)) == RealT end if RealT == Float32 @@ -221,39 +222,39 @@ isdir(outdir) && rm(outdir, recursive = true) @test_broken eltype(flux_ranocha(u_ll, u_rr, normal_direction, equations)) == RealT else - @test eltype(flux_ranocha(u_ll, u_rr, normal_direction, equations)) == RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, equations)) == RealT end - @test eltype(eltype(splitting_drikakis_tsangaris(u, normal_direction, + @test eltype(eltype(@inferred splitting_drikakis_tsangaris(u, normal_direction, equations))) == RealT - @test eltype(eltype(splitting_vanleer_haenel(u, normal_direction, equations))) == + @test eltype(eltype(@inferred splitting_vanleer_haenel(u, normal_direction, equations))) == RealT - @test eltype(eltype(splitting_lax_friedrichs(u, normal_direction, equations))) == + @test eltype(eltype(@inferred splitting_lax_friedrichs(u, normal_direction, equations))) == RealT - @test eltype(max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == + @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == RealT - @test eltype(min_max_speed_naive(u_ll, u_rr, normal_direction, equations)) == + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, equations)) == RealT - @test eltype(min_max_speed_davis(u_ll, u_rr, normal_direction, equations)) == + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction, equations)) == RealT - @test eltype(min_max_speed_einfeldt(u_ll, u_rr, normal_direction, equations)) == + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, normal_direction, equations)) == RealT for orientation in orientations - @test eltype(flux(u, orientation, equations)) == RealT - @test eltype(flux_shima_etal(u_ll, u_rr, orientation, equations)) == RealT - @test eltype(flux_kennedy_gruber(u_ll, u_rr, orientation, equations)) == + @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(flux_lmars(u_ll, u_rr, orientation, equations)) == RealT - @test eltype(flux_hllc(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(flux_chandrashekar(u_ll, u_rr, orientation, equations)) == + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, equations)) == RealT end if RealT == Float32 @@ -261,37 +262,37 @@ isdir(outdir) && rm(outdir, recursive = true) @test_broken eltype(flux_ranocha(u_ll, u_rr, orientation, equations)) == RealT else - @test eltype(flux_ranocha(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == RealT end - @test eltype(eltype(splitting_steger_warming(u, orientation, equations))) == + @test eltype(eltype(@inferred splitting_steger_warming(u, orientation, equations))) == RealT - @test eltype(eltype(splitting_drikakis_tsangaris(u, orientation, equations))) == + @test eltype(eltype(@inferred splitting_drikakis_tsangaris(u, orientation, equations))) == RealT - @test eltype(eltype(splitting_vanleer_haenel(u, orientation, equations))) == + @test eltype(eltype(@inferred splitting_vanleer_haenel(u, orientation, equations))) == RealT - @test eltype(eltype(splitting_lax_friedrichs(u, orientation, equations))) == + @test eltype(eltype(@inferred splitting_lax_friedrichs(u, orientation, equations))) == RealT - @test eltype(max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT - @test eltype(min_max_speed_naive(u_ll, u_rr, orientation, equations)) == + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, equations)) == RealT - @test eltype(min_max_speed_davis(u_ll, u_rr, orientation, equations)) == + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, equations)) == RealT - @test eltype(min_max_speed_einfeldt(u_ll, u_rr, orientation, equations)) == + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, orientation, equations)) == RealT end - @test eltype(max_abs_speeds(u, equations)) == RealT - @test eltype(cons2prim(u, equations)) == RealT - @test eltype(prim2cons(u, equations)) == RealT - @test eltype(cons2entropy(u, equations)) == RealT - @test eltype(entropy2cons(u, equations)) == RealT - # @test eltype(cons2entropy_guermond_etal(u, equations)) == RealT - @test eltype(density(u, equations)) == RealT - @test eltype(pressure(u, equations)) == RealT - @test eltype(density_pressure(u, equations)) == RealT + @test eltype(@inferred 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 + # TODO: Not implemented @test eltype(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 end end end From 77704a8f6e2ce11fc2d37ec5d4d889a49ab50893 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 27 May 2024 21:57:47 -0700 Subject: [PATCH 20/21] complete --- src/Trixi.jl | 5 +- src/equations/compressible_euler_3d.jl | 11 +- test/test_type.jl | 339 ++++++++++++++++++++----- 3 files changed, 279 insertions(+), 76 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index b001ea8846c..3c16752afa8 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -216,10 +216,11 @@ export initial_condition_eoc_test_coupled_euler_gravity, source_terms_eoc_test_coupled_euler_gravity, source_terms_eoc_test_euler export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean, - cons2entropy, entropy2cons + cons2entropy, entropy2cons, cons2entropy_guermond_etal export density, pressure, density_pressure, velocity, global_mean_vars, equilibrium_distribution, waterheight_pressure -export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, +export entropy, entropy_math, entropy_thermodynamic, entropy_guermond_etal, + energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity, enstrophy export lake_at_rest_error diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index c8476728f07..4f4d4553a8f 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -418,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 diff --git a/test/test_type.jl b/test/test_type.jl index 7e3c0519f89..bccbb50a4db 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -16,7 +16,7 @@ isdir(outdir) && rm(outdir, recursive = true) v_mean_global = (zero(RealT), zero(RealT)) c_mean_global = one(RealT) rho_mean_global = one(RealT) - equations = @inferred AcousticPerturbationEquations2D(v_mean_global, + equations = @inferred AcousticPerturbationEquations2D(v_mean_global, c_mean_global, rho_mean_global) @@ -34,11 +34,14 @@ isdir(outdir) && rm(outdir, recursive = true) 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_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 + @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 @@ -50,9 +53,9 @@ isdir(outdir) && rm(outdir, recursive = true) equations)) == RealT else @test eltype(@inferred boundary_condition_wall(u_inner, orientation, - direction, x, t, - surface_flux_function, - equations)) == RealT + direction, x, t, + surface_flux_function, + equations)) == RealT end end end @@ -64,21 +67,27 @@ isdir(outdir) && rm(outdir, recursive = true) surface_flux_function, equations)) == RealT else - @test eltype(@inferred boundary_condition_slip_wall(u_inner, normal_direction, x, t, - surface_flux_function, equations)) == + @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)) == + @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 - @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)) == + @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 - @test eltype(@inferred dissipation(u_ll, u_rr, orientation, equations)) == RealT end @test eltype(@inferred max_abs_speeds(u, equations)) == RealT @@ -103,31 +112,40 @@ isdir(outdir) && rm(outdir, recursive = true) 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_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)) == + @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_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)) == + @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_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)) == + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == RealT end if RealT == Float32 @@ -135,20 +153,28 @@ isdir(outdir) && rm(outdir, recursive = true) @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 + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT end - @test eltype(eltype(@inferred splitting_steger_warming(u, orientation, equations))) == + @test eltype(eltype(@inferred splitting_steger_warming(u, orientation, + equations))) == RealT - @test eltype(eltype(@inferred splitting_vanleer_haenel(u, orientation, equations))) == + @test eltype(eltype(@inferred splitting_vanleer_haenel(u, orientation, + equations))) == RealT - @test eltype(eltype(@inferred splitting_coirier_vanleer(u, orientation, equations))) == + @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)) == + @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 max_abs_speeds(u, equations)) == RealT @@ -172,49 +198,61 @@ isdir(outdir) && rm(outdir, recursive = true) 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)) + 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_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)) == + @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)) == + @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 + @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)) == + @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 - @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)) == + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, normal_direction, + equations)) == RealT end if RealT == Float32 @@ -222,39 +260,51 @@ isdir(outdir) && rm(outdir, recursive = true) @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 + @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))) == + 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))) == + @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)) == + @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)) == + @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)) == + @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)) == + @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)) == + @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 - @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)) == + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == RealT end if RealT == Float32 @@ -262,25 +312,34 @@ isdir(outdir) && rm(outdir, recursive = true) @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 + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT end - @test eltype(eltype(@inferred splitting_steger_warming(u, orientation, equations))) == + @test eltype(eltype(@inferred splitting_steger_warming(u, orientation, + equations))) == RealT - @test eltype(eltype(@inferred splitting_drikakis_tsangaris(u, orientation, equations))) == + @test eltype(eltype(@inferred splitting_drikakis_tsangaris(u, orientation, + equations))) == RealT - @test eltype(eltype(@inferred splitting_vanleer_haenel(u, orientation, equations))) == + @test eltype(eltype(@inferred splitting_vanleer_haenel(u, orientation, + equations))) == RealT - @test eltype(eltype(@inferred splitting_lax_friedrichs(u, orientation, equations))) == + @test eltype(eltype(@inferred splitting_lax_friedrichs(u, orientation, + equations))) == RealT - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + @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)) == + @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)) == + @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)) == + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations)) == RealT end @@ -289,10 +348,152 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred prim2cons(u, equations)) == RealT @test eltype(@inferred cons2entropy(u, equations)) == RealT @test eltype(@inferred entropy2cons(u, equations)) == RealT - # TODO: Not implemented @test eltype(cons2entropy_guermond_etal(u, equations)) == RealT + @test eltype(@inferred 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 entropy_math(cons, equations)) == RealT + @test eltype(@inferred 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 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 entropy_math(cons, equations)) == RealT + @test eltype(@inferred entropy_thermodynamic(cons, equations)) == RealT + @test eltype(@inferred energy_internal(cons, equations)) == RealT end end end From 60836bb761aab899125e064ad861fd0324c39fd2 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Fri, 31 May 2024 18:54:15 -0700 Subject: [PATCH 21/21] add Trixi --- src/Trixi.jl | 7 +++---- test/test_type.jl | 19 ++++++++++--------- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 3c16752afa8..3a882d0962c 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -181,7 +181,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal, hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, - FluxLaxFriedrichs, max_abs_speed_naive, max_abs_speeds, + FluxLaxFriedrichs, max_abs_speed_naive, FluxHLL, min_max_speed_naive, min_max_speed_davis, min_max_speed_einfeldt, FluxLMARS, FluxRotated, @@ -216,11 +216,10 @@ export initial_condition_eoc_test_coupled_euler_gravity, source_terms_eoc_test_coupled_euler_gravity, source_terms_eoc_test_euler export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean, - cons2entropy, entropy2cons, cons2entropy_guermond_etal + cons2entropy, entropy2cons export density, pressure, density_pressure, velocity, global_mean_vars, equilibrium_distribution, waterheight_pressure -export entropy, entropy_math, entropy_thermodynamic, entropy_guermond_etal, - energy_total, energy_kinetic, energy_internal, energy_magnetic, +export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity, enstrophy export lake_at_rest_error diff --git a/test/test_type.jl b/test/test_type.jl index bccbb50a4db..de02ec47110 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -90,7 +90,7 @@ isdir(outdir) && rm(outdir, recursive = true) RealT end - @test eltype(@inferred max_abs_speeds(u, 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 @@ -177,7 +177,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations)) == RealT - @test eltype(@inferred max_abs_speeds(u, 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 @@ -343,18 +343,19 @@ isdir(outdir) && rm(outdir, recursive = true) RealT end - @test eltype(@inferred max_abs_speeds(u, 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 cons2entropy_guermond_etal(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 entropy_math(cons, equations)) == RealT - @test eltype(@inferred entropy_thermodynamic(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 @@ -482,7 +483,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations)) == RealT end - @test eltype(@inferred max_abs_speeds(u, 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 @@ -491,8 +492,8 @@ isdir(outdir) && rm(outdir, recursive = true) @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 entropy_math(cons, equations)) == RealT - @test eltype(@inferred entropy_thermodynamic(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