From dbad996a3c76ccf5539b24f9dd0c873f8ad7cce2 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Fri, 17 May 2024 14:31:13 -0700 Subject: [PATCH 01/14] start main continue --- src/equations/compressible_euler_quasi_1d.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl index 6844bf9bee..85df00b24e 100644 --- a/src/equations/compressible_euler_quasi_1d.jl +++ b/src/equations/compressible_euler_quasi_1d.jl @@ -78,18 +78,19 @@ A smooth initial condition used for convergence tests in combination with """ function initial_condition_convergence_test(x, t, equations::CompressibleEulerEquationsQuasi1D) + RealT = eltype(x) c = 2 - A = 0.1 - 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 - v1 = 1.0 + v1 = 1 e = ini^2 / rho - p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2) - a = 1.5 - 0.5 * cos(x[1] * pi) + p = (equations.gamma - 1) * (e - 0.5f0 * rho * v1^2) + a = 1.5f0 - 0.5f0 * cos(x[1] * convert(RealT, pi)) return prim2cons(SVector(rho, v1, p, a), equations) end From cc38cde8506d5d907d142ce3db2bbab93e3fd9e3 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Fri, 17 May 2024 15:43:21 -0700 Subject: [PATCH 02/14] rerun CI tests --- src/equations/compressible_euler_quasi_1d.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl index 85df00b24e..c71ebb167a 100644 --- a/src/equations/compressible_euler_quasi_1d.jl +++ b/src/equations/compressible_euler_quasi_1d.jl @@ -240,17 +240,17 @@ Further details are available in the paper: # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - a_v1_avg = 0.5 * (a_ll * v1_ll + a_rr * v1_rr) - p_avg = 0.5 * (p_ll + p_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + a_v1_avg = 0.5f0 * (a_ll * v1_ll + a_rr * v1_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr) # Calculate fluxes # Ignore orientation since it is always "1" in 1D f1 = rho_mean * a_v1_avg f2 = rho_mean * a_v1_avg * v1_avg f3 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll) + 0.5f0 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll) return SVector(f1, f2, f3, zero(eltype(u_ll))) end From 92b78e82ac27997df04ec5e2b16a613898eb9a4d Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Fri, 17 May 2024 16:08:39 -0700 Subject: [PATCH 03/14] test CI tests --- src/equations/compressible_euler_quasi_1d.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl index c71ebb167a..85df00b24e 100644 --- a/src/equations/compressible_euler_quasi_1d.jl +++ b/src/equations/compressible_euler_quasi_1d.jl @@ -240,17 +240,17 @@ Further details are available in the paper: # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5f0 * (v1_ll + v1_rr) - a_v1_avg = 0.5f0 * (a_ll * v1_ll + a_rr * v1_rr) - p_avg = 0.5f0 * (p_ll + p_rr) - velocity_square_avg = 0.5f0 * (v1_ll * v1_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + a_v1_avg = 0.5 * (a_ll * v1_ll + a_rr * v1_rr) + p_avg = 0.5 * (p_ll + p_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr) # Calculate fluxes # Ignore orientation since it is always "1" in 1D f1 = rho_mean * a_v1_avg f2 = rho_mean * a_v1_avg * v1_avg f3 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5f0 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll) + 0.5 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll) return SVector(f1, f2, f3, zero(eltype(u_ll))) end From 0b9589e6d7f98e705cd7cf3b5de910f2c92ef079 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 27 May 2024 22:24:00 -0700 Subject: [PATCH 04/14] complete the first --- src/equations/compressible_euler_quasi_1d.jl | 34 ++++++++++---------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl index 85df00b24e..331f88baf2 100644 --- a/src/equations/compressible_euler_quasi_1d.jl +++ b/src/equations/compressible_euler_quasi_1d.jl @@ -81,8 +81,8 @@ function initial_condition_convergence_test(x, t, RealT = eltype(x) c = 2 A = convert(RealT, 0.1) - L = convert(RealT, 2) - f = 1 / L + L = 2 + f = 1.0f0 / L ω = 2 * convert(RealT, pi) * f ini = c + A * sin(ω * (x[1] - t)) @@ -110,18 +110,18 @@ as defined in [`initial_condition_convergence_test`](@ref). # Same settings as in `initial_condition_convergence_test`. # Derivatives calculated with ForwardDiff.jl c = 2 - A = 0.1 + A = convert(RealT, 0.1) L = 2 - f = 1 / L - ω = 2 * pi * f + f = 1.0f0 / L + ω = 2 * convert(RealT, pi) * f x1, = x ini(x1, t) = c + A * sin(ω * (x1 - t)) rho(x1, t) = ini(x1, t) - v1(x1, t) = 1.0 + v1(x1, t) = 1 e(x1, t) = ini(x1, t)^2 / rho(x1, t) - p1(x1, t) = (equations.gamma - 1) * (e(x1, t) - 0.5 * rho(x1, t) * v1(x1, t)^2) - a(x1, t) = 1.5 - 0.5 * cos(x1 * pi) + p1(x1, t) = (equations.gamma - 1) * (e(x1, t) - 0.5f0 * rho(x1, t) * v1(x1, t)^2) + a(x1, t) = 1.5f0 - 0.5f0 * cos(x1 * pi) arho(x1, t) = a(x1, t) * rho(x1, t) arhou(x1, t) = arho(x1, t) * v1(x1, t) @@ -143,7 +143,7 @@ as defined in [`initial_condition_convergence_test`](@ref). du2 = darhou_dt(x1, t) + darhouu_dx(x1, t) + a(x1, t) * dp1_dx(x1, t) du3 = daE_dt(x1, t) + dauEp_dx(x1, t) - return SVector(du1, du2, du3, 0.0) + return SVector(du1, du2, du3, 0) end # Calculate 1D flux for a single point @@ -240,17 +240,17 @@ Further details are available in the paper: # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - a_v1_avg = 0.5 * (a_ll * v1_ll + a_rr * v1_rr) - p_avg = 0.5 * (p_ll + p_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + a_v1_avg = 0.5f0 * (a_ll * v1_ll + a_rr * v1_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr) # Calculate fluxes # Ignore orientation since it is always "1" in 1D f1 = rho_mean * a_v1_avg f2 = rho_mean * a_v1_avg * v1_avg f3 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll) + 0.5f0 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll) return SVector(f1, f2, f3, zero(eltype(u_ll))) end @@ -277,13 +277,13 @@ end e_ll = a_e_ll / a_ll v1_ll = a_rho_v1_ll / a_rho_ll v_mag_ll = abs(v1_ll) - p_ll = (equations.gamma - 1) * (e_ll - 0.5 * rho_ll * v_mag_ll^2) + p_ll = (equations.gamma - 1) * (e_ll - 0.5f0 * rho_ll * v_mag_ll^2) c_ll = sqrt(equations.gamma * p_ll / rho_ll) rho_rr = a_rho_rr / a_rr e_rr = a_e_rr / a_rr v1_rr = a_rho_v1_rr / a_rho_rr v_mag_rr = abs(v1_rr) - p_rr = (equations.gamma - 1) * (e_rr - 0.5 * rho_rr * v_mag_rr^2) + p_rr = (equations.gamma - 1) * (e_rr - 0.5f0 * rho_rr * v_mag_rr^2) c_rr = sqrt(equations.gamma * p_rr / rho_rr) λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr) @@ -294,7 +294,7 @@ end rho = a_rho / a v1 = a_rho_v1 / a_rho e = a_e / a - p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2) + p = (equations.gamma - 1) * (e - 0.5f0 * rho * v1^2) c = sqrt(equations.gamma * p / rho) return (abs(v1) + c,) From 8aa7542b35353a477c2096f2cc81e68d851dd753 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Wed, 5 Jun 2024 16:38:48 -0700 Subject: [PATCH 05/14] compressible euler quasi 1d --- src/equations/compressible_euler_quasi_1d.jl | 9 ++-- test/test_type.jl | 47 ++++++++++++++++++++ 2 files changed, 51 insertions(+), 5 deletions(-) diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl index 331f88baf2..5ae1cf4228 100644 --- a/src/equations/compressible_euler_quasi_1d.jl +++ b/src/equations/compressible_euler_quasi_1d.jl @@ -109,6 +109,7 @@ as defined in [`initial_condition_convergence_test`](@ref). equations::CompressibleEulerEquationsQuasi1D) # Same settings as in `initial_condition_convergence_test`. # Derivatives calculated with ForwardDiff.jl + RealT = eltype(u) c = 2 A = convert(RealT, 0.1) L = 2 @@ -158,7 +159,7 @@ end f2 = a_rho_v1 * v1 f3 = a * v1 * (e + p) - return SVector(f1, f2, f3, zero(eltype(u))) + return SVector(f1, f2, f3, 0) end """ @@ -190,9 +191,7 @@ Further details are available in the paper: # in the arithmetic average of {p}. p_avg = p_ll + p_rr - z = zero(eltype(u_ll)) - - return SVector(z, a_ll * p_avg, z, z) + return SVector(0, a_ll * p_avg, 0, 0) end # While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that @@ -252,7 +251,7 @@ Further details are available in the paper: f3 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + 0.5f0 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll) - return SVector(f1, f2, f3, zero(eltype(u_ll))) + return SVector(f1, f2, f3, 0) end # While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that diff --git a/test/test_type.jl b/test/test_type.jl index de02ec4711..b9828aa911 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -497,6 +497,53 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred energy_internal(cons, equations)) == RealT end end + + @timed_testset "Compressible Euler Quasi 1D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleEulerEquationsQuasi1D(RealT(1.4)) + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = SVector(one(RealT), one(RealT), one(RealT), one(RealT)) + orientation = 1 + normal_direction = normal_ll = normal_rr = SVector(one(RealT)) + + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_nonconservative_chan_etal(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_nonconservative_chan_etal(u_ll, u_rr, + normal_direction, + equations)) == + RealT + @test eltype(@inferred flux_nonconservative_chan_etal(u_ll, u_rr, normal_ll, + normal_rr, equations)) == + RealT + if RealT == Float32 + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(flux_chan_etal(u_ll, u_rr, orientation, equations)) == + RealT + else + @test eltype(@inferred flux_chan_etal(u_ll, u_rr, orientation, equations)) == + RealT + end + + @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + # TODO: There is a bug in the `entropy` function + # @test eltype(@inferred entropy(u, equations)) == RealT + @test eltype(@inferred cons2entropy(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 end # module From 6f2b008381413f2978dd30f33e02cf198414aec0 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Wed, 5 Jun 2024 22:13:16 -0700 Subject: [PATCH 06/14] compressible euler multicomponent 1d --- .../compressible_euler_multicomponent_1d.jl | 133 +++++++++--------- 1 file changed, 68 insertions(+), 65 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl index 6338e04c3e..e964284f4c 100644 --- a/src/equations/compressible_euler_multicomponent_1d.jl +++ b/src/equations/compressible_euler_multicomponent_1d.jl @@ -123,14 +123,15 @@ A smooth initial condition used for convergence tests in combination with """ function initial_condition_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations1D) + RealT = eltype(x) c = 2 - A = 0.1 + A = convert(RealT, 0.1) L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f ini = c + A * sin(omega * (x[1] - t)) - v1 = 1.0 + v1 = 1 rho = ini @@ -159,20 +160,21 @@ Source terms used for convergence tests in combination with @inline function source_terms_convergence_test(u, x, t, equations::CompressibleEulerMulticomponentEquations1D) # Same settings as in `initial_condition` + RealT = eltype(u) c = 2 - A = 0.1 + A = convert(RealT, 0.1) L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f gamma = totalgamma(u, equations) x1, = x si, co = sincos((t - x1) * omega) - tmp = (-((4 * si * A - 4c) + 1) * (gamma - 1) * co * A * omega) / 2 + tmp = (-((4 * si * A - 4 * c) + 1) * (gamma - 1) * co * A * omega) / 2 # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1 - du_rho = SVector{ncomponents(equations), real(equations)}(0.0 + du_rho = SVector{ncomponents(equations), real(equations)}(0 for i in eachcomponent(equations)) du1 = tmp @@ -194,24 +196,25 @@ A for multicomponent adapted weak blast wave adapted to multicomponent and taken function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerMulticomponentEquations1D) # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) - inicenter = SVector(0.0) + RealT = eltype(x) + inicenter = SVector(0) x_norm = x[1] - inicenter[1] r = abs(x_norm) - cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm) + cos_phi = x_norm > 0 ? 1 : -1 - prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5 ? + prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5f0 ? 2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * - 1.0 : + 1 : 2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * - 1.1691 + convert(RealT, 1.1691) for i in eachcomponent(equations)) - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - p = r > 0.5 ? 1.0 : 1.245 + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) prim_other = SVector{2, real(equations)}(v1, p) @@ -227,7 +230,7 @@ end v1 = rho_v1 / rho gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * v1^2) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * v1^2) f_rho = densities(u, v1, equations) f1 = rho_v1 * v1 + p @@ -255,7 +258,7 @@ Entropy conserving two-point flux by rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 2], u_rr[i + 2]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 2] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 2] + u_rr[i + 2]) for i in eachcomponent(equations)) @@ -269,13 +272,13 @@ Entropy conserving two-point flux by # extract velocities v1_ll = rho_v1_ll / rho_ll v1_rr = rho_v1_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) - v1_square = 0.5 * (v1_ll^2 + v1_rr^2) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v1_square = 0.5f0 * (v1_ll^2 + v1_rr^2) v_sum = v1_avg - enth = zero(v_sum) - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) + enth = 0 + help1_ll = 0 + help1_rr = 0 for i in eachcomponent(equations) enth += rhok_avg[i] * gas_constants[i] @@ -283,14 +286,14 @@ Entropy conserving two-point flux by help1_rr += u_rr[i + 2] * cv[i] end - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2)) / help1_rr - T = 0.5 * (1.0 / T_ll + 1.0 / T_rr) - T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) + T_ll = (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2)) / help1_ll + T_rr = (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2)) / help1_rr + T = 0.5f0 * (1 / T_ll + 1 / T_rr) + T_log = ln_mean(1 / T_ll, 1 / T_rr) # Calculate fluxes depending on orientation - help1 = zero(T_ll) - help2 = zero(T_rr) + help1 = 0 + help2 = 0 f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg for i in eachcomponent(equations)) @@ -299,7 +302,7 @@ Entropy conserving two-point flux by help2 += f_rho[i] end f1 = (help2) * v1_avg + enth / T - f2 = (help1) / T_log - 0.5 * (v1_square) * (help2) + v1_avg * f1 + f2 = (help1) / T_log - 0.5f0 * (v1_square) * (help2) + v1_avg * f1 f_other = SVector{2, real(equations)}(f1, f2) @@ -330,7 +333,7 @@ See also rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 2], u_rr[i + 2]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 2] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 2] + u_rr[i + 2]) for i in eachcomponent(equations)) @@ -339,25 +342,25 @@ See also rho_rr = density(u_rr, equations) # Calculating gamma - gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + gamma = totalgamma(0.5f0 * (u_ll + u_rr), equations) inv_gamma_minus_one = 1 / (gamma - 1) # extract velocities v1_ll = rho_v1_ll / rho_ll v1_rr = rho_v1_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr) # density flux f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg for i in eachcomponent(equations)) # helpful variables - f_rho_sum = zero(v1_ll) - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) - enth_ll = zero(v1_ll) - enth_rr = zero(v1_rr) + f_rho_sum = 0 + help1_ll = 0 + help1_rr = 0 + enth_ll = 0 + enth_rr = 0 for i in eachcomponent(equations) enth_ll += u_ll[i + 2] * gas_constants[i] enth_rr += u_rr[i + 2] * gas_constants[i] @@ -367,17 +370,17 @@ See also end # temperature and pressure - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2)) / help1_rr + T_ll = (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2)) / help1_ll + T_rr = (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2)) / help1_rr p_ll = T_ll * enth_ll p_rr = T_rr * enth_rr - p_avg = 0.5 * (p_ll + p_rr) + p_avg = 0.5f0 * (p_ll + p_rr) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) # momentum and energy flux f1 = f_rho_sum * v1_avg + p_avg f2 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (p_ll * v1_rr + p_rr * v1_ll) + 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) f_other = SVector{2, real(equations)}(f1, f2) return vcat(f_other, f_rho) @@ -398,8 +401,8 @@ end v_ll = rho_v1_ll / rho_ll v_rr = rho_v1_rr / rho_rr - p_ll = (gamma_ll - 1) * (rho_e_ll - 1 / 2 * rho_ll * v_ll^2) - p_rr = (gamma_rr - 1) * (rho_e_rr - 1 / 2 * rho_rr * v_rr^2) + p_ll = (gamma_ll - 1) * (rho_e_ll - 0.5f0 * rho_ll * v_ll^2) + p_rr = (gamma_rr - 1) * (rho_e_rr - 0.5f0 * rho_rr * v_rr^2) c_ll = sqrt(gamma_ll * p_ll / rho_ll) c_rr = sqrt(gamma_rr * p_rr / rho_rr) @@ -414,7 +417,7 @@ end v1 = rho_v1 / rho gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 1 / 2 * rho * (v1^2)) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2)) c = sqrt(gamma * p / rho) return (abs(v1) + c,) @@ -431,7 +434,7 @@ end v1 = rho_v1 / rho gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2)) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2)) prim_other = SVector{2, real(equations)}(v1, p) return vcat(prim_other, prim_rho) @@ -451,7 +454,7 @@ end rho_v1 = rho * v1 - rho_e = p / (gamma - 1) + 0.5 * (rho_v1 * v1) + rho_e = p / (gamma - 1) + 0.5f0 * (rho_v1 * v1) cons_other = SVector{2, RealT}(rho_v1, rho_e) @@ -466,8 +469,8 @@ end rho = density(u, equations) - help1 = zero(rho) - gas_constant = zero(rho) + help1 = 0 + gas_constant = 0 for i in eachcomponent(equations) help1 += u[i + 2] * cv[i] gas_constant += gas_constants[i] * (u[i + 2] / rho) @@ -477,10 +480,10 @@ end v_square = v1^2 gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * v_square) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * v_square) s = log(p) - gamma * log(rho) - log(gas_constant) rho_p = rho / p - T = (rho_e - 0.5 * rho * v_square) / (help1) + T = (rho_e - 0.5f0 * rho * v_square) / (help1) entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] * (1 - log(T)) + @@ -507,14 +510,14 @@ end (-cv[i] * log(-w[2]) - cp[i] + w[i + 2] - - 0.5 * w[1]^2 / + 0.5f0 * w[1]^2 / w[2])) for i in eachcomponent(equations)) - rho = zero(cons_rho[1]) - help1 = zero(cons_rho[1]) - help2 = zero(cons_rho[1]) - p = zero(cons_rho[1]) + rho = 0 + help1 = 0 + help2 = 0 + p = 0 for i in eachcomponent(equations) rho += cons_rho[i] help1 += cons_rho[i] * cv[i] * gammas[i] @@ -523,7 +526,7 @@ end end u1 = rho * v1 gamma = help1 / help2 - u2 = p / (gamma - 1) + 0.5 * rho * v1^2 + u2 = p / (gamma - 1) + 0.5f0 * rho * v1^2 cons_other = SVector{2, real(equations)}(u1, u2) return vcat(cons_other, cons_rho) end @@ -534,7 +537,7 @@ end rho = density(u, equations) T = temperature(u, equations) - total_entropy = zero(u[1]) + total_entropy = 0 for i in eachcomponent(equations) total_entropy -= u[i + 2] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 2])) end @@ -548,7 +551,7 @@ end rho_v1, rho_e = u rho = density(u, equations) - help1 = zero(rho) + help1 = 0 for i in eachcomponent(equations) help1 += u[i + 2] * cv[i] @@ -556,7 +559,7 @@ end v1 = rho_v1 / rho v_square = v1^2 - T = (rho_e - 0.5 * rho * v_square) / help1 + T = (rho_e - 0.5f0 * rho * v_square) / help1 return T end @@ -570,8 +573,8 @@ partial density fractions as well as the partial specific heats at constant volu @inline function totalgamma(u, equations::CompressibleEulerMulticomponentEquations1D) @unpack cv, gammas = equations - help1 = zero(u[1]) - help2 = zero(u[1]) + help1 = 0 + help2 = 0 for i in eachcomponent(equations) help1 += u[i + 2] * cv[i] * gammas[i] @@ -587,13 +590,13 @@ end rho = density(u, equations) gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * (rho_v1^2) / rho) + p = (gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2) / rho) return p end @inline function density(u, equations::CompressibleEulerMulticomponentEquations1D) - rho = zero(u[1]) + rho = 0 for i in eachcomponent(equations) rho += u[i + 2] From 0cecb0efbc20e0d5eae7a9850c43748513b90103 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Wed, 5 Jun 2024 23:01:17 -0700 Subject: [PATCH 07/14] compressible euler multicomponent 2d --- .../compressible_euler_multicomponent_2d.jl | 173 +++++++++--------- 1 file changed, 88 insertions(+), 85 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 60fce222f2..5e9c027ca6 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -127,15 +127,16 @@ A smooth initial condition used for convergence tests in combination with """ function initial_condition_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations2D) + RealT = eltype(x) c = 2 - A = 0.1 + A = convert(RealT, 0.1) L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f ini = c + A * sin(omega * (x[1] + x[2] - t)) - v1 = 1.0 - v2 = 1.0 + v1 = 1 + v2 = 1 rho = ini @@ -165,11 +166,12 @@ Source terms used for convergence tests in combination with @inline function source_terms_convergence_test(u, x, t, equations::CompressibleEulerMulticomponentEquations2D) # Same settings as in `initial_condition` + RealT = eltype(u) c = 2 - A = 0.1 + A = convert(RealT, 0.1) L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f gamma = totalgamma(u, equations) @@ -191,7 +193,7 @@ Source terms used for convergence tests in combination with du1 = tmp5 du2 = tmp5 - du3 = 2 * ((tmp6 - 1.0) * tmp3 + tmp6 * gamma) * tmp1 + du3 = 2 * ((tmp6 - 1) * tmp3 + tmp6 * gamma) * tmp1 du_other = SVector{3, real(equations)}(du1, du2, du3) @@ -210,27 +212,28 @@ function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerMulticomponentEquations2D) # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # Set up polar coordinates - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) phi = atan(y_norm, x_norm) sin_phi, cos_phi = sincos(phi) - prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5 ? + prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5f0 ? 2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * - 1.0 : + 1 : 2^(i - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * - 1.1691 + convert(RealT, 1.1691) for i in eachcomponent(equations)) - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi - p = r > 0.5 ? 1.0 : 1.245 + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) prim_other = SVector{3, real(equations)}(v1, v2, p) @@ -247,7 +250,7 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2 + v2^2)) if orientation == 1 f_rho = densities(u, v1, equations) @@ -277,7 +280,7 @@ end v2 = rho_v2 / rho v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2 + v2^2)) f_rho = densities(u, v_normal, equations) f1 = rho_v1 * v_normal + p * normal_direction[1] @@ -306,7 +309,7 @@ Adaption of the entropy conserving two-point flux by rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], u_rr[i + 3]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 3] + u_rr[i + 3]) for i in eachcomponent(equations)) @@ -319,15 +322,15 @@ Adaption of the entropy conserving two-point flux by v2_ll = rho_v2_ll / rho_ll v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v1_square = 0.5 * (v1_ll^2 + v1_rr^2) - v2_square = 0.5 * (v2_ll^2 + v2_rr^2) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v1_square = 0.5f0 * (v1_ll^2 + v1_rr^2) + v2_square = 0.5f0 * (v2_ll^2 + v2_rr^2) v_sum = v1_avg + v2_avg - enth = zero(v_sum) - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) + enth = 0 + help1_ll = 0 + help1_rr = 0 for i in eachcomponent(equations) enth += rhok_avg[i] * gas_constants[i] @@ -335,14 +338,14 @@ Adaption of the entropy conserving two-point flux by help1_rr += u_rr[i + 3] * cv[i] end - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr - T = 0.5 * (1.0 / T_ll + 1.0 / T_rr) - T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) + T_ll = (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll + T_rr = (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr + T = 0.5f0 * (1 / T_ll + 1 / T_rr) + T_log = ln_mean(1 / T_ll, 1 / T_rr) # Calculate fluxes depending on orientation - help1 = zero(T_ll) - help2 = zero(T_rr) + help1 = 0 + help2 = 0 if orientation == 1 f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg for i in eachcomponent(equations)) @@ -352,7 +355,7 @@ Adaption of the entropy conserving two-point flux by end f1 = (help2) * v1_avg + enth / T f2 = (help2) * v2_avg - f3 = (help1) / T_log - 0.5 * (v1_square + v2_square) * (help2) + v1_avg * f1 + + f3 = (help1) / T_log - 0.5f0 * (v1_square + v2_square) * (help2) + v1_avg * f1 + v2_avg * f2 else f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v2_avg @@ -363,7 +366,7 @@ Adaption of the entropy conserving two-point flux by end f1 = (help2) * v1_avg f2 = (help2) * v2_avg + enth / T - f3 = (help1) / T_log - 0.5 * (v1_square + v2_square) * (help2) + v1_avg * f1 + + f3 = (help1) / T_log - 0.5f0 * (v1_square + v2_square) * (help2) + v1_avg * f1 + v2_avg * f2 end f_other = SVector{3, real(equations)}(f1, f2, f3) @@ -395,7 +398,7 @@ See also rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], u_rr[i + 3]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 3] + u_rr[i + 3]) for i in eachcomponent(equations)) @@ -404,23 +407,23 @@ See also rho_rr = density(u_rr, equations) # Calculating gamma - gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + gamma = totalgamma(0.5f0 * (u_ll + u_rr), equations) inv_gamma_minus_one = 1 / (gamma - 1) # extract velocities v1_ll = rho_v1_ll / rho_ll v1_rr = rho_v1_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) v2_ll = rho_v2_ll / rho_ll v2_rr = rho_v2_rr / rho_rr - v2_avg = 0.5 * (v2_ll + v2_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) # helpful variables - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) - enth_ll = zero(v1_ll) - enth_rr = zero(v1_rr) + help1_ll = 0 + help1_rr = 0 + enth_ll = 0 + enth_rr = 0 for i in eachcomponent(equations) enth_ll += u_ll[i + 3] * gas_constants[i] enth_rr += u_rr[i + 3] * gas_constants[i] @@ -429,14 +432,14 @@ See also end # temperature and pressure - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr + T_ll = (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll + T_rr = (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr p_ll = T_ll * enth_ll p_rr = T_rr * enth_rr - p_avg = 0.5 * (p_ll + p_rr) + p_avg = 0.5f0 * (p_ll + p_rr) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - f_rho_sum = zero(T_rr) + f_rho_sum = 0 if orientation == 1 f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg for i in eachcomponent(equations)) @@ -446,7 +449,7 @@ See also f1 = f_rho_sum * v1_avg + p_avg f2 = f_rho_sum * v2_avg f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (p_ll * v1_rr + p_rr * v1_ll) + 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) else f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v2_avg for i in eachcomponent(equations)) @@ -456,7 +459,7 @@ See also f1 = f_rho_sum * v1_avg f2 = f_rho_sum * v2_avg + p_avg f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (p_ll * v2_rr + p_rr * v2_ll) + 0.5f0 * (p_ll * v2_rr + p_rr * v2_ll) end # momentum and energy flux @@ -474,7 +477,7 @@ end rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], u_rr[i + 3]) for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 3] + u_rr[i + 3]) for i in eachcomponent(equations)) @@ -483,25 +486,25 @@ end rho_rr = density(u_rr, equations) # Calculating gamma - gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + gamma = totalgamma(0.5f0 * (u_ll + u_rr), equations) inv_gamma_minus_one = 1 / (gamma - 1) # extract velocities v1_ll = rho_v1_ll / rho_ll v1_rr = rho_v1_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) v2_ll = rho_v2_ll / rho_ll v2_rr = rho_v2_rr / rho_rr - v2_avg = 0.5 * (v2_ll + v2_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] # helpful variables - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) - enth_ll = zero(v1_ll) - enth_rr = zero(v1_rr) + help1_ll = 0 + help1_rr = 0 + enth_ll = 0 + enth_rr = 0 for i in eachcomponent(equations) enth_ll += u_ll[i + 3] * gas_constants[i] enth_rr += u_rr[i + 3] * gas_constants[i] @@ -510,15 +513,15 @@ end end # temperature and pressure - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr + T_ll = (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll + T_rr = (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr p_ll = T_ll * enth_ll p_rr = T_rr * enth_rr - p_avg = 0.5 * (p_ll + p_rr) + p_avg = 0.5f0 * (p_ll + p_rr) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - f_rho_sum = zero(T_rr) - f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5 * + f_rho_sum = 0 + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) for i in eachcomponent(equations)) for i in eachcomponent(equations) @@ -527,7 +530,7 @@ end f1 = f_rho_sum * v1_avg + p_avg * normal_direction[1] f2 = f_rho_sum * v2_avg + p_avg * normal_direction[2] f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) + 0.5f0 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) # momentum and energy flux f_other = SVector(f1, f2, f3) @@ -557,9 +560,9 @@ end end # Compute the sound speeds on the left and right - p_ll = (gamma_ll - 1) * (rho_e_ll - 1 / 2 * (rho_v1_ll^2 + rho_v2_ll^2) / rho_ll) + p_ll = (gamma_ll - 1) * (rho_e_ll - 0.5f0 * (rho_v1_ll^2 + rho_v2_ll^2) / rho_ll) c_ll = sqrt(gamma_ll * p_ll / rho_ll) - p_rr = (gamma_rr - 1) * (rho_e_rr - 1 / 2 * (rho_v1_rr^2 + rho_v2_rr^2) / rho_rr) + p_rr = (gamma_rr - 1) * (rho_e_rr - 0.5f0 * (rho_v1_rr^2 + rho_v2_rr^2) / rho_rr) c_rr = sqrt(gamma_rr * p_rr / rho_rr) λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) @@ -574,7 +577,7 @@ end v2 = rho_v2 / rho gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 1 / 2 * rho * (v1^2 + v2^2)) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2 + v2^2)) c = sqrt(gamma * p / rho) return (abs(v1) + c, abs(v2) + c) @@ -635,7 +638,7 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2 + v2^2)) prim_other = SVector{3, real(equations)}(v1, v2, p) return vcat(prim_other, prim_rho) @@ -649,8 +652,8 @@ end rho = density(u, equations) # Multicomponent stuff - help1 = zero(rho) - gas_constant = zero(rho) + help1 = 0 + gas_constant = 0 for i in eachcomponent(equations) help1 += u[i + 3] * cv[i] gas_constant += gas_constants[i] * (u[i + 3] / rho) @@ -661,10 +664,10 @@ end v_square = v1^2 + v2^2 gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * v_square) + p = (gamma - 1) * (rho_e - 0.5f0 * rho * v_square) s = log(p) - gamma * log(rho) - log(gas_constant) rho_p = rho / p - T = (rho_e - 0.5 * rho * v_square) / (help1) + T = (rho_e - 0.5f0 * rho * v_square) / (help1) entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] * (1 - log(T)) + @@ -698,10 +701,10 @@ end 1) for i in eachcomponent(equations)) - rho = zero(cons_rho[1]) - help1 = zero(cons_rho[1]) - help2 = zero(cons_rho[1]) - p = zero(cons_rho[1]) + rho = 0 + help1 = 0 + help2 = 0 + p = 0 for i in eachcomponent(equations) rho += cons_rho[i] help1 += cons_rho[i] * cv[i] * gammas[i] @@ -711,7 +714,7 @@ end u1 = rho * v1 u2 = rho * v2 gamma = help1 / help2 - u3 = p / (gamma - 1) + 0.5 * rho * v_squared + u3 = p / (gamma - 1) + 0.5f0 * rho * v_squared cons_other = SVector{3, real(equations)}(u1, u2, u3) return vcat(cons_other, cons_rho) end @@ -728,7 +731,7 @@ end rho_v1 = rho * v1 rho_v2 = rho * v2 - rho_e = p / (gamma - 1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2) + rho_e = p / (gamma - 1) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2) cons_other = SVector{3, real(equations)}(rho_v1, rho_v2, rho_e) @@ -740,7 +743,7 @@ end rho = density(u, equations) T = temperature(u, equations) - total_entropy = zero(u[1]) + total_entropy = 0 for i in eachcomponent(equations) total_entropy -= u[i + 3] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 3])) end @@ -754,7 +757,7 @@ end rho_v1, rho_v2, rho_e = u rho = density(u, equations) - help1 = zero(rho) + help1 = 0 for i in eachcomponent(equations) help1 += u[i + 3] * cv[i] @@ -763,7 +766,7 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v_square = v1^2 + v2^2 - T = (rho_e - 0.5 * rho * v_square) / help1 + T = (rho_e - 0.5f0 * rho * v_square) / help1 return T end @@ -777,8 +780,8 @@ partial density fractions as well as the partial specific heats at constant volu @inline function totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D) @unpack cv, gammas = equations - help1 = zero(u[1]) - help2 = zero(u[1]) + help1 = 0 + help2 = 0 for i in eachcomponent(equations) help1 += u[i + 3] * cv[i] * gammas[i] @@ -794,13 +797,13 @@ end rho = density(u, equations) gamma = totalgamma(u, equations) - rho_times_p = (gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2)) + rho_times_p = (gamma - 1) * (rho * rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2)) return rho_times_p end @inline function density(u, equations::CompressibleEulerMulticomponentEquations2D) - rho = zero(u[1]) + rho = 0 for i in eachcomponent(equations) rho += u[i + 3] From 2249e591066776f8b922cce86035db2910190165 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 6 Jun 2024 00:51:59 -0700 Subject: [PATCH 08/14] complete unit tests --- test/test_type.jl | 167 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 144 insertions(+), 23 deletions(-) diff --git a/test/test_type.jl b/test/test_type.jl index b9828aa911..bdd37f5f55 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -47,10 +47,13 @@ isdir(outdir) && rm(outdir, recursive = true) for direction in directions if RealT == Float32 # check `surface_flux_function` (test broken) - @test_broken eltype(boundary_condition_wall(u_inner, orientation, - direction, x, t, - surface_flux_function, - equations)) == RealT + @test_broken eltype(@inferred boundary_condition_wall(u_inner, + orientation, + direction, x, + t, + surface_flux_function, + equations)) == + RealT else @test eltype(@inferred boundary_condition_wall(u_inner, orientation, direction, x, t, @@ -62,10 +65,12 @@ isdir(outdir) && rm(outdir, recursive = true) if RealT == Float32 # check `surface_flux_function` (test broken) - @test_broken eltype(boundary_condition_slip_wall(u_inner, normal_direction, - x, t, - surface_flux_function, - equations)) == RealT + @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner, + normal_direction, + x, t, + surface_flux_function, + equations)) == + RealT else @test eltype(@inferred boundary_condition_slip_wall(u_inner, normal_direction, x, t, @@ -141,7 +146,8 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred flux_hllc(u_ll, u_rr, orientation, equations)) == RealT if RealT == Float32 # check `ln_mean` (test broken) - @test_broken eltype(flux_chandrashekar(u_ll, u_rr, orientation, equations)) == + @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == RealT else @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, @@ -150,7 +156,8 @@ isdir(outdir) && rm(outdir, recursive = true) end if RealT == Float32 # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(flux_ranocha(u_ll, u_rr, orientation, equations)) == + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, + equations)) == RealT else @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == @@ -247,8 +254,9 @@ isdir(outdir) && rm(outdir, recursive = true) RealT if RealT == Float32 # check `ln_mean` (test broken) - @test_broken eltype(flux_chandrashekar(u_ll, u_rr, normal_direction, - equations)) == + @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, + normal_direction, + equations)) == RealT else @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, normal_direction, @@ -257,7 +265,8 @@ isdir(outdir) && rm(outdir, recursive = true) end if RealT == Float32 # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(flux_ranocha(u_ll, u_rr, normal_direction, equations)) == + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, + equations)) == RealT else @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, @@ -299,8 +308,9 @@ isdir(outdir) && rm(outdir, recursive = true) RealT if RealT == Float32 # check `ln_mean` (test broken) - @test_broken eltype(flux_chandrashekar(u_ll, u_rr, orientation, - equations)) == + @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, + orientation, + equations)) == RealT else @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, @@ -309,7 +319,8 @@ isdir(outdir) && rm(outdir, recursive = true) end if RealT == Float32 # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(flux_ranocha(u_ll, u_rr, orientation, equations)) == + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, + equations)) == RealT else @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == @@ -357,6 +368,7 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred Trixi.entropy_math(cons, equations)) == RealT @test eltype(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT @test eltype(@inferred energy_internal(cons, equations)) == RealT + # TODO: test `gradient_conservative`, not necessary but good to have end end @@ -416,15 +428,17 @@ isdir(outdir) && rm(outdir, recursive = true) RealT if RealT == Float32 # check `ln_mean` (test broken) - @test_broken eltype(flux_chandrashekar(u_ll, u_rr, normal_direction, - equations)) == RealT + @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, + normal_direction, + equations)) == RealT else @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, normal_direction, equations)) == RealT end if RealT == Float32 # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(flux_ranocha(u_ll, u_rr, normal_direction, equations)) == + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, + equations)) == RealT else @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, @@ -454,15 +468,17 @@ isdir(outdir) && rm(outdir, recursive = true) RealT if RealT == Float32 # check `ln_mean` (test broken) - @test_broken eltype(flux_chandrashekar(u_ll, u_rr, orientation, - equations)) == RealT + @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, + orientation, + equations)) == RealT else @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, equations)) == RealT end if RealT == Float32 # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(flux_ranocha(u_ll, u_rr, orientation, equations)) == + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, + equations)) == RealT else @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == @@ -498,9 +514,113 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Compressible Euler Multicomponent 1D" begin + for RealT in (Float32, Float64) + gammas = (RealT(1.4), RealT(1.4)) + gas_constants = (RealT(0.4), RealT(0.4)) + equations = @inferred CompressibleEulerMulticomponentEquations1D(gammas = gammas, + gas_constants = gas_constants) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = SVector(one(RealT), one(RealT), one(RealT), one(RealT)) + orientation = 1 + + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + # Note: Pass with check for `ln_mean` + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, equations)) == + RealT + # Note: Pass with check for `ln_mean` and `inv_ln_mean` + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT + + @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test eltype(@inferred entropy2cons(u, equations)) == RealT + @test eltype(@inferred Trixi.total_entropy(u, equations)) == RealT + @test eltype(@inferred Trixi.temperature(u, equations)) == RealT + @test eltype(@inferred Trixi.totalgamma(u, equations)) == RealT + @test eltype(@inferred density(u, equations)) == RealT + @test eltype(@inferred pressure(u, equations)) == RealT + end + end + + @timed_testset "Compressible Euler Multicomponent 2D" begin + for RealT in (Float32, Float64) + gammas = (RealT(1.4), RealT(1.4)) + gas_constants = (RealT(0.4), RealT(0.4)) + equations = @inferred CompressibleEulerMulticomponentEquations2D(gammas = gammas, + gas_constants = gas_constants) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = SVector(one(RealT), one(RealT), one(RealT), one(RealT), + one(RealT)) + orientations = [1, 2] + normal_direction = SVector(one(RealT), zero(RealT)) + + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, + equations)) == RealT + else + @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, + equations)) == RealT + end + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + # pass with check for `ln_mean` + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == + RealT + # pass with check for `ln_mean` and `inv_ln_mean` + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT + + @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == + RealT + end + + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test eltype(@inferred entropy2cons(u, equations)) == RealT + @test eltype(@inferred Trixi.total_entropy(u, equations)) == RealT + @test eltype(@inferred Trixi.temperature(u, equations)) == RealT + @test eltype(@inferred Trixi.totalgamma(u, equations)) == RealT + @test eltype(@inferred density(u, equations)) == RealT + @test eltype(@inferred density_pressure(u, equations)) == RealT + end + end + @timed_testset "Compressible Euler Quasi 1D" begin for RealT in (Float32, Float64) equations = @inferred CompressibleEulerEquationsQuasi1D(RealT(1.4)) + x = SVector(zero(RealT)) t = zero(RealT) u = u_ll = u_rr = SVector(one(RealT), one(RealT), one(RealT), one(RealT)) @@ -524,7 +644,8 @@ isdir(outdir) && rm(outdir, recursive = true) RealT if RealT == Float32 # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(flux_chan_etal(u_ll, u_rr, orientation, equations)) == + @test_broken eltype(@inferred flux_chan_etal(u_ll, u_rr, orientation, + equations)) == RealT else @test eltype(@inferred flux_chan_etal(u_ll, u_rr, orientation, equations)) == From b1c167bdc5a349f13c4ef9d62922b31056ea340b Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Sat, 8 Jun 2024 21:26:22 -0700 Subject: [PATCH 09/14] fix comments --- test/test_type.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_type.jl b/test/test_type.jl index bdd37f5f55..7507eacfbf 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -535,10 +535,10 @@ isdir(outdir) && rm(outdir, recursive = true) RealT @test eltype(@inferred flux(u, orientation, equations)) == RealT - # Note: Pass with check for `ln_mean` + # pass with check for `ln_mean` @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, equations)) == RealT - # Note: Pass with check for `ln_mean` and `inv_ln_mean` + # pass with check for `ln_mean` and `inv_ln_mean` @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == RealT @@ -581,7 +581,7 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred flux(u, normal_direction, equations)) == RealT if RealT == Float32 - # check `ln_mean` (test broken) + # check `ln_mean` and `inv_ln_mean` (test broken) @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, equations)) == RealT else From 930815cf7e34e223e7991256e104d59a077955b5 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 10 Jun 2024 21:58:09 -0700 Subject: [PATCH 10/14] complete --- .../compressible_euler_multicomponent_1d.jl | 69 ++++++++------- .../compressible_euler_multicomponent_2d.jl | 83 ++++++++++--------- test/test_type.jl | 50 +++++++---- 3 files changed, 120 insertions(+), 82 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl index e964284f4c..d4f6b421f7 100644 --- a/src/equations/compressible_euler_multicomponent_1d.jl +++ b/src/equations/compressible_euler_multicomponent_1d.jl @@ -145,7 +145,7 @@ function initial_condition_convergence_test(x, t, prim1 = rho * v1 prim2 = rho^2 - prim_other = SVector{2, real(equations)}(prim1, prim2) + prim_other = SVector(prim1, prim2) return vcat(prim_other, prim_rho) end @@ -180,7 +180,7 @@ Source terms used for convergence tests in combination with du1 = tmp du2 = tmp - du_other = SVector{2, real(equations)}(du1, du2) + du_other = SVector(du1, du2) return vcat(du_other, du_rho) end @@ -216,7 +216,7 @@ function initial_condition_weak_blast_wave(x, t, v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) - prim_other = SVector{2, real(equations)}(v1, p) + prim_other = SVector(v1, p) return prim2cons(vcat(prim_other, prim_rho), equations) end @@ -236,7 +236,7 @@ end f1 = rho_v1 * v1 + p f2 = (rho_e + p) * v1 - f_other = SVector{2, real(equations)}(f1, f2) + f_other = SVector(f1, f2) return vcat(f_other, f_rho) end @@ -276,9 +276,10 @@ Entropy conserving two-point flux by v1_square = 0.5f0 * (v1_ll^2 + v1_rr^2) v_sum = v1_avg - enth = 0 - help1_ll = 0 - help1_rr = 0 + RealT = eltype(u_ll) + enth = zero(RealT) + help1_ll = zero(RealT) + help1_rr = zero(RealT) for i in eachcomponent(equations) enth += rhok_avg[i] * gas_constants[i] @@ -292,8 +293,8 @@ Entropy conserving two-point flux by T_log = ln_mean(1 / T_ll, 1 / T_rr) # Calculate fluxes depending on orientation - help1 = 0 - help2 = 0 + help1 = zero(RealT) + help2 = zero(RealT) f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg for i in eachcomponent(equations)) @@ -304,7 +305,7 @@ Entropy conserving two-point flux by f1 = (help2) * v1_avg + enth / T f2 = (help1) / T_log - 0.5f0 * (v1_square) * (help2) + v1_avg * f1 - f_other = SVector{2, real(equations)}(f1, f2) + f_other = SVector(f1, f2) return vcat(f_other, f_rho) end @@ -356,11 +357,12 @@ See also for i in eachcomponent(equations)) # helpful variables - f_rho_sum = 0 - help1_ll = 0 - help1_rr = 0 - enth_ll = 0 - enth_rr = 0 + RealT = eltype(u_ll) + f_rho_sum = zero(RealT) + help1_ll = zero(RealT) + help1_rr = zero(RealT) + enth_ll = zero(RealT) + enth_rr = zero(RealT) for i in eachcomponent(equations) enth_ll += u_ll[i + 2] * gas_constants[i] enth_rr += u_rr[i + 2] * gas_constants[i] @@ -381,7 +383,7 @@ See also f1 = f_rho_sum * v1_avg + p_avg f2 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) - f_other = SVector{2, real(equations)}(f1, f2) + f_other = SVector(f1, f2) return vcat(f_other, f_rho) end @@ -435,7 +437,7 @@ end gamma = totalgamma(u, equations) p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2)) - prim_other = SVector{2, real(equations)}(v1, p) + prim_other = SVector(v1, p) return vcat(prim_other, prim_rho) end @@ -469,8 +471,9 @@ end rho = density(u, equations) - help1 = 0 - gas_constant = 0 + RealT = eltype(u) + help1 = zero(RealT) + gas_constant = zero(RealT) for i in eachcomponent(equations) help1 += u[i + 2] * cv[i] gas_constant += gas_constants[i] * (u[i + 2] / rho) @@ -495,7 +498,7 @@ end w1 = gas_constant * v1 * rho_p w2 = gas_constant * (-rho_p) - entrop_other = SVector{2, real(equations)}(w1, w2) + entrop_other = SVector(w1, w2) return vcat(entrop_other, entrop_rho) end @@ -514,10 +517,11 @@ end w[2])) for i in eachcomponent(equations)) - rho = 0 - help1 = 0 - help2 = 0 - p = 0 + RealT = eltype(w) + rho = zero(RealT) + help1 = zero(RealT) + help2 = zero(RealT) + p = zero(RealT) for i in eachcomponent(equations) rho += cons_rho[i] help1 += cons_rho[i] * cv[i] * gammas[i] @@ -527,7 +531,7 @@ end u1 = rho * v1 gamma = help1 / help2 u2 = p / (gamma - 1) + 0.5f0 * rho * v1^2 - cons_other = SVector{2, real(equations)}(u1, u2) + cons_other = SVector(u1, u2) return vcat(cons_other, cons_rho) end @@ -537,7 +541,8 @@ end rho = density(u, equations) T = temperature(u, equations) - total_entropy = 0 + RealT = eltype(u) + total_entropy = zero(RealT) for i in eachcomponent(equations) total_entropy -= u[i + 2] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 2])) end @@ -551,7 +556,9 @@ end rho_v1, rho_e = u rho = density(u, equations) - help1 = 0 + + RealT = eltype(u) + help1 = zero(RealT) for i in eachcomponent(equations) help1 += u[i + 2] * cv[i] @@ -573,8 +580,9 @@ partial density fractions as well as the partial specific heats at constant volu @inline function totalgamma(u, equations::CompressibleEulerMulticomponentEquations1D) @unpack cv, gammas = equations - help1 = 0 - help2 = 0 + RealT = eltype(u) + help1 = zero(RealT) + help2 = zero(RealT) for i in eachcomponent(equations) help1 += u[i + 2] * cv[i] * gammas[i] @@ -596,7 +604,8 @@ end end @inline function density(u, equations::CompressibleEulerMulticomponentEquations1D) - rho = 0 + RealT = eltype(u) + rho = zero(RealT) for i in eachcomponent(equations) rho += u[i + 2] diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 5e9c027ca6..3473f88733 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -151,7 +151,7 @@ function initial_condition_convergence_test(x, t, prim2 = rho * v2 prim3 = rho^2 - prim_other = SVector{3, real(equations)}(prim1, prim2, prim3) + prim_other = SVector(prim1, prim2, prim3) return vcat(prim_other, prim_rho) end @@ -195,7 +195,7 @@ Source terms used for convergence tests in combination with du2 = tmp5 du3 = 2 * ((tmp6 - 1) * tmp3 + tmp6 * gamma) * tmp1 - du_other = SVector{3, real(equations)}(du1, du2, du3) + du_other = SVector(du1, du2, du3) return vcat(du_other, du_rho) end @@ -235,7 +235,7 @@ function initial_condition_weak_blast_wave(x, t, v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) - prim_other = SVector{3, real(equations)}(v1, v2, p) + prim_other = SVector(v1, v2, p) return prim2cons(vcat(prim_other, prim_rho), equations) end @@ -264,7 +264,7 @@ end f3 = (rho_e + p) * v2 end - f_other = SVector{3, real(equations)}(f1, f2, f3) + f_other = SVector(f1, f2, f3) return vcat(f_other, f_rho) end @@ -287,7 +287,7 @@ end f2 = rho_v2 * v_normal + p * normal_direction[2] f3 = (rho_e + p) * v_normal - f_other = SVector{3, real(equations)}(f1, f2, f3) + f_other = SVector(f1, f2, f3) return vcat(f_other, f_rho) end @@ -328,9 +328,10 @@ Adaption of the entropy conserving two-point flux by v2_square = 0.5f0 * (v2_ll^2 + v2_rr^2) v_sum = v1_avg + v2_avg - enth = 0 - help1_ll = 0 - help1_rr = 0 + RealT = eltype(u_ll) + enth = zero(RealT) + help1_ll = zero(RealT) + help1_rr = zero(RealT) for i in eachcomponent(equations) enth += rhok_avg[i] * gas_constants[i] @@ -344,8 +345,8 @@ Adaption of the entropy conserving two-point flux by T_log = ln_mean(1 / T_ll, 1 / T_rr) # Calculate fluxes depending on orientation - help1 = 0 - help2 = 0 + help1 = zero(RealT) + help2 = zero(RealT) if orientation == 1 f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg for i in eachcomponent(equations)) @@ -369,7 +370,7 @@ Adaption of the entropy conserving two-point flux by f3 = (help1) / T_log - 0.5f0 * (v1_square + v2_square) * (help2) + v1_avg * f1 + v2_avg * f2 end - f_other = SVector{3, real(equations)}(f1, f2, f3) + f_other = SVector(f1, f2, f3) return vcat(f_other, f_rho) end @@ -420,10 +421,11 @@ See also velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) # helpful variables - help1_ll = 0 - help1_rr = 0 - enth_ll = 0 - enth_rr = 0 + RealT = eltype(u_ll) + help1_ll = zero(RealT) + help1_rr = zero(RealT) + enth_ll = zero(RealT) + enth_rr = zero(RealT) for i in eachcomponent(equations) enth_ll += u_ll[i + 3] * gas_constants[i] enth_rr += u_rr[i + 3] * gas_constants[i] @@ -439,7 +441,7 @@ See also p_avg = 0.5f0 * (p_ll + p_rr) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - f_rho_sum = 0 + f_rho_sum = zero(RealT) if orientation == 1 f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg for i in eachcomponent(equations)) @@ -463,7 +465,7 @@ See also end # momentum and energy flux - f_other = SVector{3, real(equations)}(f1, f2, f3) + f_other = SVector(f1, f2, f3) return vcat(f_other, f_rho) end @@ -501,10 +503,11 @@ end v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] # helpful variables - help1_ll = 0 - help1_rr = 0 - enth_ll = 0 - enth_rr = 0 + RealT = eltype(u_ll) + help1_ll = zero(RealT) + help1_rr = zero(RealT) + enth_ll = zero(RealT) + enth_rr = zero(RealT) for i in eachcomponent(equations) enth_ll += u_ll[i + 3] * gas_constants[i] enth_rr += u_rr[i + 3] * gas_constants[i] @@ -520,7 +523,7 @@ end p_avg = 0.5f0 * (p_ll + p_rr) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - f_rho_sum = 0 + f_rho_sum = zero(RealT) f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) for i in eachcomponent(equations)) @@ -639,7 +642,7 @@ end v2 = rho_v2 / rho gamma = totalgamma(u, equations) p = (gamma - 1) * (rho_e - 0.5f0 * rho * (v1^2 + v2^2)) - prim_other = SVector{3, real(equations)}(v1, v2, p) + prim_other = SVector(v1, v2, p) return vcat(prim_other, prim_rho) end @@ -652,8 +655,9 @@ end rho = density(u, equations) # Multicomponent stuff - help1 = 0 - gas_constant = 0 + RealT = eltype(u) + help1 = zero(RealT) + gas_constant = zero(RealT) for i in eachcomponent(equations) help1 += u[i + 3] * cv[i] gas_constant += gas_constants[i] * (u[i + 3] / rho) @@ -680,7 +684,7 @@ end w2 = gas_constant * v2 * rho_p w3 = gas_constant * (-rho_p) - entrop_other = SVector{3, real(equations)}(w1, w2, w3) + entrop_other = SVector(w1, w2, w3) return vcat(entrop_other, entrop_rho) end @@ -701,10 +705,11 @@ end 1) for i in eachcomponent(equations)) - rho = 0 - help1 = 0 - help2 = 0 - p = 0 + RealT = eltype(w) + rho = zero(RealT) + help1 = zero(RealT) + help2 = zero(RealT) + p = zero(RealT) for i in eachcomponent(equations) rho += cons_rho[i] help1 += cons_rho[i] * cv[i] * gammas[i] @@ -715,7 +720,7 @@ end u2 = rho * v2 gamma = help1 / help2 u3 = p / (gamma - 1) + 0.5f0 * rho * v_squared - cons_other = SVector{3, real(equations)}(u1, u2, u3) + cons_other = SVector(u1, u2, u3) return vcat(cons_other, cons_rho) end @@ -733,7 +738,7 @@ end rho_v2 = rho * v2 rho_e = p / (gamma - 1) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2) - cons_other = SVector{3, real(equations)}(rho_v1, rho_v2, rho_e) + cons_other = SVector(rho_v1, rho_v2, rho_e) return vcat(cons_other, cons_rho) end @@ -743,7 +748,8 @@ end rho = density(u, equations) T = temperature(u, equations) - total_entropy = 0 + RealT = eltype(u) + total_entropy = zero(RealT) for i in eachcomponent(equations) total_entropy -= u[i + 3] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 3])) end @@ -757,7 +763,8 @@ end rho_v1, rho_v2, rho_e = u rho = density(u, equations) - help1 = 0 + RealT = eltype(u) + help1 = zero(RealT) for i in eachcomponent(equations) help1 += u[i + 3] * cv[i] @@ -780,8 +787,9 @@ partial density fractions as well as the partial specific heats at constant volu @inline function totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D) @unpack cv, gammas = equations - help1 = 0 - help2 = 0 + RealT = eltype(u) + help1 = zero(RealT) + help2 = zero(RealT) for i in eachcomponent(equations) help1 += u[i + 3] * cv[i] * gammas[i] @@ -803,7 +811,8 @@ end end @inline function density(u, equations::CompressibleEulerMulticomponentEquations2D) - rho = 0 + RealT = eltype(u) + rho = zero(RealT) for i in eachcomponent(equations) rho += u[i + 3] diff --git a/test/test_type.jl b/test/test_type.jl index 7507eacfbf..41d3b62893 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -535,12 +535,23 @@ isdir(outdir) && rm(outdir, recursive = true) RealT @test eltype(@inferred flux(u, orientation, equations)) == RealT - # pass with check for `ln_mean` - @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, equations)) == - RealT - # pass with check for `ln_mean` and `inv_ln_mean` - @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == - RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == RealT + else + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == RealT + end + if RealT == Float32 + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, + equations)) == + RealT + else + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT + end @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @@ -591,13 +602,23 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations @test eltype(@inferred flux(u, orientation, equations)) == RealT - # pass with check for `ln_mean` - @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, - equations)) == - RealT - # pass with check for `ln_mean` and `inv_ln_mean` - @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == - RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, + orientation, + equations)) == RealT + else + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == RealT + end + if RealT == Float32 + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, + equations)) == RealT + else + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT + end @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == @@ -657,8 +678,7 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT @test eltype(@inferred cons2prim(u, equations)) == RealT @test eltype(@inferred prim2cons(u, equations)) == RealT - # TODO: There is a bug in the `entropy` function - # @test eltype(@inferred entropy(u, equations)) == RealT + @test eltype(@inferred entropy(u, equations)) == RealT @test eltype(@inferred cons2entropy(u, equations)) == RealT @test eltype(@inferred density(u, equations)) == RealT @test eltype(@inferred pressure(u, equations)) == RealT From b6bd36e291cda5c65bfaecac4f159d9ae87a7f86 Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Thu, 13 Jun 2024 20:58:10 -1000 Subject: [PATCH 11/14] apply strict type check Co-authored-by: Hendrik Ranocha --- test/test_type.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/test/test_type.jl b/test/test_type.jl index 41d3b62893..4617fa9abd 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -630,11 +630,11 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred prim2cons(u, equations)) == RealT @test eltype(@inferred cons2entropy(u, equations)) == RealT @test eltype(@inferred entropy2cons(u, equations)) == RealT - @test eltype(@inferred Trixi.total_entropy(u, equations)) == RealT - @test eltype(@inferred Trixi.temperature(u, equations)) == RealT - @test eltype(@inferred Trixi.totalgamma(u, equations)) == RealT - @test eltype(@inferred density(u, equations)) == RealT - @test eltype(@inferred density_pressure(u, equations)) == RealT + @test typeof(@inferred Trixi.total_entropy(u, equations)) == RealT + @test typeof(@inferred Trixi.temperature(u, equations)) == RealT + @test typeof(@inferred Trixi.totalgamma(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT end end @@ -678,11 +678,11 @@ isdir(outdir) && rm(outdir, recursive = true) @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 entropy(u, equations)) == RealT + @test typeof(@inferred entropy(u, equations)) == RealT @test eltype(@inferred cons2entropy(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 typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT end end end From c74b5a479738a0b56c091b187644f83f2b4ac87d Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 13 Jun 2024 21:21:46 -1000 Subject: [PATCH 12/14] cover all --- test/test_type.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/test/test_type.jl b/test/test_type.jl index 4617fa9abd..748efba62c 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -553,18 +553,18 @@ isdir(outdir) && rm(outdir, recursive = true) RealT end - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT @test eltype(@inferred cons2prim(u, equations)) == RealT @test eltype(@inferred prim2cons(u, equations)) == RealT @test eltype(@inferred cons2entropy(u, equations)) == RealT @test eltype(@inferred entropy2cons(u, equations)) == RealT - @test eltype(@inferred Trixi.total_entropy(u, equations)) == RealT - @test eltype(@inferred Trixi.temperature(u, equations)) == RealT - @test eltype(@inferred Trixi.totalgamma(u, equations)) == RealT - @test eltype(@inferred density(u, equations)) == RealT - @test eltype(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred Trixi.total_entropy(u, equations)) == RealT + @test typeof(@inferred Trixi.temperature(u, equations)) == RealT + @test typeof(@inferred Trixi.totalgamma(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT end end @@ -620,7 +620,7 @@ isdir(outdir) && rm(outdir, recursive = true) RealT end - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT end @@ -673,13 +673,13 @@ isdir(outdir) && rm(outdir, recursive = true) RealT end - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT @test eltype(@inferred cons2prim(u, equations)) == RealT @test eltype(@inferred prim2cons(u, equations)) == RealT - @test typeof(@inferred entropy(u, equations)) == RealT @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred entropy(u, equations)) == RealT @test typeof(@inferred density(u, equations)) == RealT @test typeof(@inferred pressure(u, equations)) == RealT @test typeof(@inferred density_pressure(u, equations)) == RealT From 307c2ec159f814589a0ca8f62d40185bec8168fd Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 13 Jun 2024 21:47:16 -1000 Subject: [PATCH 13/14] cover before --- test/test_type.jl | 54 +++++++++++++++++++++++------------------------ 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/test/test_type.jl b/test/test_type.jl index 748efba62c..7bea104bb8 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -80,7 +80,7 @@ isdir(outdir) && rm(outdir, recursive = true) end @test eltype(@inferred flux(u, normal_direction, equations)) == RealT - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == RealT @test eltype(@inferred dissipation(u_ll, u_rr, normal_direction, equations)) == @@ -88,7 +88,7 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations @test eltype(@inferred flux(u, orientation, equations)) == RealT - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @test eltype(@inferred dissipation(u_ll, u_rr, orientation, equations)) == @@ -174,7 +174,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations))) == RealT - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @@ -189,11 +189,11 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred prim2cons(u, equations)) == RealT @test eltype(@inferred cons2entropy(u, equations)) == RealT @test eltype(@inferred entropy2cons(u, equations)) == RealT - @test eltype(@inferred density(u, equations)) == RealT - @test eltype(@inferred pressure(u, equations)) == RealT - @test eltype(@inferred density_pressure(u, equations)) == RealT - @test eltype(@inferred entropy(cons, equations)) == RealT - @test eltype(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT end end @@ -282,7 +282,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations))) == RealT - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == RealT @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, @@ -340,7 +340,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations))) == RealT - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, @@ -359,15 +359,15 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred prim2cons(u, equations)) == RealT @test eltype(@inferred cons2entropy(u, equations)) == RealT @test eltype(@inferred entropy2cons(u, equations)) == RealT - @test eltype(@inferred Trixi.entropy_guermond_etal(u, equations)) == RealT @test eltype(@inferred Trixi.cons2entropy_guermond_etal(u, equations)) == RealT - @test eltype(@inferred density(u, equations)) == RealT - @test eltype(@inferred pressure(u, equations)) == RealT - @test eltype(@inferred density_pressure(u, equations)) == RealT - @test eltype(@inferred entropy(cons, equations)) == RealT - @test eltype(@inferred Trixi.entropy_math(cons, equations)) == RealT - @test eltype(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT - @test eltype(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_guermond_etal(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_math(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT # TODO: test `gradient_conservative`, not necessary but good to have end @@ -445,7 +445,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations)) == RealT end - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == RealT @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, @@ -489,7 +489,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations))) == RealT - @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @@ -504,13 +504,13 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred prim2cons(u, equations)) == RealT @test eltype(@inferred cons2entropy(u, equations)) == RealT @test eltype(@inferred entropy2cons(u, equations)) == RealT - @test eltype(@inferred density(u, equations)) == RealT - @test eltype(@inferred pressure(u, equations)) == RealT - @test eltype(@inferred density_pressure(u, equations)) == RealT - @test eltype(@inferred entropy(cons, equations)) == RealT - @test eltype(@inferred Trixi.entropy_math(cons, equations)) == RealT - @test eltype(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT - @test eltype(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_math(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT end end From 71b7a647803e7266a4063a0d8424941355d4a667 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 13 Jun 2024 21:52:10 -1000 Subject: [PATCH 14/14] add me --- AUTHORS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/AUTHORS.md b/AUTHORS.md index 8e3afcf867..54d6321633 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -43,3 +43,4 @@ are listed in alphabetical order: * Michael Schlottke-Lakemper * Toskan Theine * Andrew Winters +* Huiyu Xie