Skip to content

Commit

Permalink
HLL Davis wave speeds & Cartesian Flux Winters for polytropic Euler (#…
Browse files Browse the repository at this point in the history
…1733)

* hll_davis for polytropic euler

* flux winters for tree and tests

* fmt

* typo

* remove copied comment

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Andrew Winters <[email protected]>

* orientation_or_normal_direction

* typo

---------

Co-authored-by: Andrew Winters <[email protected]>
  • Loading branch information
DanielDoehring and andrewwinters5000 authored Nov 17, 2023
1 parent 8bdd0cd commit c8e1f99
Show file tree
Hide file tree
Showing 6 changed files with 272 additions and 2 deletions.
63 changes: 63 additions & 0 deletions examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@

using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the polytropic Euler equations

gamma = 1.4
kappa = 0.5 # Scaling factor for the pressure.
equations = PolytropicEulerEquations2D(gamma, kappa)

initial_condition = initial_condition_convergence_test

volume_flux = flux_winters_etal
solver = DGSEM(polydeg = 3, surface_flux = FluxHLL(min_max_speed_davis),
volume_integral = VolumeIntegralFluxDifferencing(volume_flux))

coordinates_min = (0.0, 0.0)
coordinates_max = (1.0, 1.0)

# Create a uniformly refined mesh with periodic boundaries
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 2,
periodicity = true,
n_cells_max = 30_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
source_terms = source_terms_convergence_test)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 0.1)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
extra_analysis_errors = (:l2_error_primitive,
:linf_error_primitive))

alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(interval = 100,
save_initial_solution = true,
save_final_solution = true,
solution_variables = cons2prim)

stepsize_callback = StepsizeCallback(cfl = 0.1)

callbacks = CallbackSet(summary_callback,
analysis_callback, alive_callback,
save_solution,
stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
save_everystep = false, callback = callbacks);
summary_callback() # print the timer summary
102 changes: 100 additions & 2 deletions src/equations/polytropic_euler_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ function initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquat
return prim2cons(SVector(rho, v1, v2), equations)
end

# Calculate 1D flux for a single point in the normal direction
# Calculate 2D flux for a single point in the normal direction
# Note, this directional vector is not normalized
@inline function flux(u, normal_direction::AbstractVector,
equations::PolytropicEulerEquations2D)
Expand All @@ -135,8 +135,28 @@ end
return SVector(f1, f2, f3)
end

# Calculate 2D flux for a single point
@inline function flux(u, orientation::Integer, equations::PolytropicEulerEquations2D)
_, v1, v2 = cons2prim(u, equations)
p = pressure(u, equations)

rho_v1 = u[2]
rho_v2 = u[3]

if orientation == 1
f1 = rho_v1
f2 = rho_v1 * v1 + p
f3 = rho_v1 * v2
else
f1 = rho_v2
f2 = rho_v2 * v1
f3 = rho_v2 * v2 + p
end
return SVector(f1, f2, f3)
end

"""
flux_winters_etal(u_ll, u_rr, normal_direction,
flux_winters_etal(u_ll, u_rr, orientation_or_normal_direction,
equations::PolytropicEulerEquations2D)
Entropy conserving two-point flux for isothermal or polytropic gases.
Expand Down Expand Up @@ -178,6 +198,37 @@ For details see Section 3.2 of the following reference
return SVector(f1, f2, f3)
end

@inline function flux_winters_etal(u_ll, u_rr, orientation::Integer,
equations::PolytropicEulerEquations2D)
# Unpack left and right state
rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
p_ll = equations.kappa * rho_ll^equations.gamma
p_rr = equations.kappa * rho_rr^equations.gamma

# Compute the necessary mean values
if equations.gamma == 1.0 # 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)

if orientation == 1 # x-direction
f1 = rho_mean * 0.5 * (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)
f2 = f1 * v1_avg
f3 = f1 * v2_avg + p_avg
end

return SVector(f1, f2, f3)
end

@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
equations::PolytropicEulerEquations2D)
rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
Expand All @@ -196,6 +247,53 @@ end
return lambda_min, lambda_max
end

# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
equations::PolytropicEulerEquations2D)
rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
# Pressure for polytropic Euler
p_ll = equations.kappa * rho_ll^equations.gamma
p_rr = equations.kappa * rho_rr^equations.gamma

c_ll = sqrt(equations.gamma * p_ll / rho_ll)
c_rr = sqrt(equations.gamma * p_rr / rho_rr)

if orientation == 1 # x-direction
λ_min = min(v1_ll - c_ll, v1_rr - c_rr)
λ_max = max(v1_ll + c_ll, v1_rr + c_rr)
else # y-direction
λ_min = min(v2_ll - c_ll, v2_rr - c_rr)
λ_max = max(v2_ll + c_ll, v2_rr + c_rr)
end

return λ_min, λ_max
end

# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
@inline function min_max_speed_davis(u_ll, u_rr, normal_direction::AbstractVector,
equations::PolytropicEulerEquations2D)
rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
# Pressure for polytropic Euler
p_ll = equations.kappa * rho_ll^equations.gamma
p_rr = equations.kappa * rho_rr^equations.gamma

norm_ = norm(normal_direction)

c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_
c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_

v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]

# The v_normals are already scaled by the norm
λ_min = min(v_normal_ll - c_ll, v_normal_rr - c_rr)
λ_max = max(v_normal_ll + c_ll, v_normal_rr + c_rr)

return λ_min, λ_max
end

@inline function max_abs_speeds(u, equations::PolytropicEulerEquations2D)
rho, v1, v2 = cons2prim(u, equations)
c = sqrt(equations.gamma * equations.kappa * rho^(equations.gamma - 1))
Expand Down
23 changes: 23 additions & 0 deletions test/test_structured_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -598,6 +598,29 @@ end
end
end

@trixi_testset "elixir_eulerpolytropic_convergence.jl: HLL(Davis)" begin
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"),
solver=DGSEM(polydeg = 3,
surface_flux = FluxHLL(min_max_speed_davis),
volume_integral = VolumeIntegralFluxDifferencing(volume_flux)),
l2=[
0.0016689832177644243, 0.0025920263793104445,
0.003281074494629298,
],
linf=[
0.01099488320190023, 0.013309526619350365,
0.02008032661117909,
])
# Ensure that we do not have excessive memory allocations
# (e.g., from type instabilities)
let
t = sol.t[end]
u_ode = sol.u[end]
du_ode = similar(u_ode)
@test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
end
end

@trixi_testset "elixir_eulerpolytropic_ec.jl" begin
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_ec.jl"),
l2=[
Expand Down
35 changes: 35 additions & 0 deletions test/test_tree_2d_eulerpolytropic.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
module TestExamples2DEulerMulticomponent

using Test
using Trixi

include("test_trixi.jl")

EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem")

@testset "Polytropic Euler" begin
#! format: noindent

@trixi_testset "elixir_eulerpolytropic_convergence.jl" begin
@test_trixi_include(joinpath(EXAMPLES_DIR,
"elixir_eulerpolytropic_convergence.jl"),
l2=[
0.0016689832177626373, 0.0025920263793094526,
0.003281074494626679,
],
linf=[
0.010994883201896677, 0.013309526619350365,
0.02008032661117376,
])
# Ensure that we do not have excessive memory allocations
# (e.g., from type instabilities)
let
t = sol.t[end]
u_ode = sol.u[end]
du_ode = similar(u_ode)
@test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
end
end
end

end # module
3 changes: 3 additions & 0 deletions test/test_tree_2d_part2.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,9 @@ isdir(outdir) && rm(outdir, recursive = true)
# Compressible Euler Multicomponent
include("test_tree_2d_eulermulti.jl")

# Compressible Polytropic Euler
include("test_tree_2d_eulerpolytropic.jl")

# Compressible Euler coupled with acoustic perturbation equations
include("test_tree_2d_euleracoustics.jl")

Expand Down
48 changes: 48 additions & 0 deletions test/test_unit.jl
Original file line number Diff line number Diff line change
Expand Up @@ -782,6 +782,54 @@ end
end
end

@timed_testset "Consistency check for HLL flux with Davis wave speed estimates: Polytropic CEE" begin
flux_hll = FluxHLL(min_max_speed_davis)

gamma = 1.4
kappa = 0.5 # Scaling factor for the pressure.
equations = PolytropicEulerEquations2D(gamma, kappa)
u = SVector(1.1, -0.5, 2.34)

orientations = [1, 2]
for orientation in orientations
@test flux_hll(u, u, orientation, equations) flux(u, orientation, equations)
end

normal_directions = [SVector(1.0, 0.0),
SVector(0.0, 1.0),
SVector(0.5, -0.5),
SVector(-1.2, 0.3)]

for normal_direction in normal_directions
@test flux_hll(u, u, normal_direction, equations)
flux(u, normal_direction, equations)
end
end

@timed_testset "Consistency check for Winters flux: Polytropic CEE" begin
for gamma in [1.4, 1.0, 5 / 3]
kappa = 0.5 # Scaling factor for the pressure.
equations = PolytropicEulerEquations2D(gamma, kappa)
u = SVector(1.1, -0.5, 2.34)

orientations = [1, 2]
for orientation in orientations
@test flux_winters_etal(u, u, orientation, equations)
flux(u, orientation, equations)
end

normal_directions = [SVector(1.0, 0.0),
SVector(0.0, 1.0),
SVector(0.5, -0.5),
SVector(-1.2, 0.3)]

for normal_direction in normal_directions
@test flux_winters_etal(u, u, normal_direction, equations)
flux(u, normal_direction, equations)
end
end
end

@timed_testset "Consistency check for HLL flux with Davis wave speed estimates: LEE" begin
flux_hll = FluxHLL(min_max_speed_davis)

Expand Down

0 comments on commit c8e1f99

Please sign in to comment.