Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add numerical support of other real types (polytropic_euler) #2015

Merged
merged 8 commits into from
Aug 16, 2024
Merged

Add numerical support of other real types (polytropic_euler) #2015

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
82f2294
start
huiyuxie Jul 19, 2024
9f52433
Merge branch 'main' into main-polytropic
huiyuxie Jul 24, 2024
0cf0589
Merge branch 'main' into main-polytropic
huiyuxie Jul 29, 2024
d151122
fix flux tests
huiyuxie Jul 29, 2024
c6fecc4
keep conflicts
huiyuxie Jul 29, 2024
21b2012
complete equations
huiyuxie Jul 29, 2024
d900e3c
Merge branch 'main' into main-polytropic
huiyuxie Aug 15, 2024
eff115a
complete unit test
huiyuxie Aug 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 27 additions & 25 deletions src/equations/polytropic_euler_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ in combination with [`source_terms_convergence_test`](@ref).
function initial_condition_convergence_test(x, t, equations::PolytropicEulerEquations2D)
# manufactured initial condition from Winters (2019) [0.1007/s10543-019-00789-w]
# domain must be set to [0, 1] x [0, 1]
h = 8 + cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t)
h = 8 + cospi(2 * x[1]) * sinpi(2 * x[2]) * cospi(2 * t)

return SVector(h, h / 2, 3 * h / 2)
end
Expand All @@ -78,19 +78,20 @@ Source terms used for convergence tests in combination with
rho, v1, v2 = cons2prim(u, equations)

# Residual from Winters (2019) [0.1007/s10543-019-00789-w] eq. (5.2).
h = 8 + cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t)
h_t = -2 * pi * cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * sin(2 * pi * t)
h_x = -2 * pi * sin(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t)
h_y = 2 * pi * cos(2 * pi * x[1]) * cos(2 * pi * x[2]) * cos(2 * pi * t)
RealT = eltype(u)
h = 8 + cospi(2 * x[1]) * sinpi(2 * x[2]) * cospi(2 * t)
h_t = -2 * convert(RealT, pi) * cospi(2 * x[1]) * sinpi(2 * x[2]) * sinpi(2 * t)
h_x = -2 * convert(RealT, pi) * sinpi(2 * x[1]) * sinpi(2 * x[2]) * cospi(2 * t)
h_y = 2 * convert(RealT, pi) * cospi(2 * x[1]) * cospi(2 * x[2]) * cospi(2 * t)

rho_x = h_x
rho_y = h_y

b = equations.kappa * equations.gamma * h^(equations.gamma - 1)

r_1 = h_t + h_x / 2 + 3 / 2 * h_y
r_2 = h_t / 2 + h_x / 4 + b * rho_x + 3 / 4 * h_y
r_3 = 3 / 2 * h_t + 3 / 4 * h_x + 9 / 4 * h_y + b * rho_y
r_1 = h_t + h_x / 2 + 3 * h_y / 2
r_2 = h_t / 2 + h_x / 4 + b * rho_x + 3 * h_y / 4
r_3 = 3 * h_t / 2 + 3 * h_x / 4 + 9 * h_y / 4 + b * rho_y

return SVector(r_1, r_2, r_3)
end
Expand All @@ -113,9 +114,10 @@ function initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquat
phi = atan(y_norm, x_norm)

# 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)
RealT = eltype(x)
rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691)
v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi)
v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin(phi)

return prim2cons(SVector(rho, v1, v2), equations)
end
Expand Down Expand Up @@ -181,17 +183,17 @@ For details see Section 3.2 of the following reference
v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]

# Compute the necessary mean values
if equations.gamma == 1.0 # isothermal gas
if equations.gamma == 1 # isothermal gas
rho_mean = ln_mean(rho_ll, rho_rr)
else # equations.gamma > 1 # polytropic gas
rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma)
end
v1_avg = 0.5 * (v1_ll + v1_rr)
v2_avg = 0.5 * (v2_ll + v2_rr)
p_avg = 0.5 * (p_ll + p_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)

# 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]

Expand All @@ -207,21 +209,21 @@ end
p_rr = equations.kappa * rho_rr^equations.gamma

# Compute the necessary mean values
if equations.gamma == 1.0 # isothermal gas
if equations.gamma == 1 # isothermal gas
rho_mean = ln_mean(rho_ll, rho_rr)
else # equations.gamma > 1 # polytropic gas
rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma)
end
v1_avg = 0.5 * (v1_ll + v1_rr)
v2_avg = 0.5 * (v2_ll + v2_rr)
p_avg = 0.5 * (p_ll + p_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)

if orientation == 1 # x-direction
f1 = rho_mean * 0.5 * (v1_ll + v1_rr)
f1 = rho_mean * 0.5f0 * (v1_ll + v1_rr)
f2 = f1 * v1_avg + p_avg
f3 = f1 * v2_avg
else # y-direction
f1 = rho_mean * 0.5 * (v2_ll + v2_rr)
f1 = rho_mean * 0.5f0 * (v2_ll + v2_rr)
f2 = f1 * v1_avg
f3 = f1 * v2_avg + p_avg
end
Expand Down Expand Up @@ -360,14 +362,14 @@ end
v_square = v1^2 + v2^2
p = pressure(u, equations)
# Form of the internal energy depends on gas type
if equations.gamma == 1.0 # isothermal gas
if equations.gamma == 1 # isothermal gas
internal_energy = equations.kappa * log(rho)
else # equations.gamma > 1 # polytropic gas
internal_energy = equations.kappa * rho^(equations.gamma - 1) /
(equations.gamma - 1.0)
(equations.gamma - 1)
end

w1 = internal_energy + p / rho - 0.5 * v_square
w1 = internal_energy + p / rho - 0.5f0 * v_square
w2 = v1
w3 = v2

Expand Down
74 changes: 74 additions & 0 deletions test/test_type.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1926,6 +1926,80 @@ isdir(outdir) && rm(outdir, recursive = true)
end
end

@timed_testset "Polytropic Euler 2D" begin
for RealT in (Float32, Float64)
equations1 = @inferred PolytropicEulerEquations2D(RealT(1),
RealT(1)) # equations.gamma == 1
equations2 = @inferred PolytropicEulerEquations2D(RealT(1.4), RealT(0.5)) # equations.gamma > 1

for equations in (equations1, equations2)
x = SVector(zero(RealT), zero(RealT))
t = zero(RealT)
u = u_ll = u_rr = SVector(one(RealT), one(RealT), one(RealT))
orientations = [1, 2]
directions = [1, 2, 3, 4]
normal_direction = SVector(one(RealT), zero(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 initial_condition_weak_blast_wave(x, t, equations)) ==
RealT

@test eltype(@inferred flux(u, normal_direction, equations)) == RealT
if RealT == Float32
# check `ln_mean` and `stolarsky_mean` (test broken)
@test_broken eltype(@inferred flux_winters_etal(u_ll, u_rr,
normal_direction,
equations)) ==
RealT
else
@test eltype(@inferred flux_winters_etal(u_ll, u_rr, normal_direction,
equations)) ==
RealT
end
@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 typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction,
equations)) ==
RealT

for orientation in orientations
@test eltype(@inferred flux(u, orientation, equations)) == RealT
if RealT == Float32
# check `ln_mean` and `stolarsky_mean` (test broken)
@test_broken eltype(@inferred flux_winters_etal(u_ll, u_rr,
orientation,
equations)) ==
RealT
else
@test eltype(@inferred flux_winters_etal(u_ll, u_rr, orientation,
equations)) ==
RealT
end
@test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation,
equations)) ==
RealT
@test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation,
equations)) ==
RealT
end

@test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
@test eltype(@inferred cons2prim(u, equations)) == RealT
@test eltype(@inferred prim2cons(u, equations)) == RealT
@test eltype(@inferred cons2entropy(u, equations)) == RealT
@test typeof(@inferred density(u, equations)) == RealT
@test typeof(@inferred pressure(u, equations)) == RealT
end
end
end

@timed_testset "Shallow Water 1D" begin
for RealT in (Float32, Float64)
equations = @inferred ShallowWaterEquations1D(gravity_constant = RealT(9.81))
Expand Down
Loading