From c8e1f99dc52666cee2d5d8000d881e7f3dac60e0 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 17 Nov 2023 08:36:37 +0100 Subject: [PATCH 1/3] HLL Davis wave speeds & Cartesian Flux Winters for polytropic Euler (#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 * orientation_or_normal_direction * typo --------- Co-authored-by: Andrew Winters --- .../elixir_eulerpolytropic_convergence.jl | 63 +++++++++++ src/equations/polytropic_euler_2d.jl | 102 +++++++++++++++++- test/test_structured_2d.jl | 23 ++++ test/test_tree_2d_eulerpolytropic.jl | 35 ++++++ test/test_tree_2d_part2.jl | 3 + test/test_unit.jl | 48 +++++++++ 6 files changed, 272 insertions(+), 2 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl create mode 100644 test/test_tree_2d_eulerpolytropic.jl diff --git a/examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl b/examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl new file mode 100644 index 00000000000..95ef38eb701 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl @@ -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 diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index d4902bbafb2..f5d2f7b0bad 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -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) @@ -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. @@ -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) @@ -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)) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index dd2248e10b2..17ba001d5ec 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -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=[ diff --git a/test/test_tree_2d_eulerpolytropic.jl b/test/test_tree_2d_eulerpolytropic.jl new file mode 100644 index 00000000000..545cf7274ff --- /dev/null +++ b/test/test_tree_2d_eulerpolytropic.jl @@ -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 diff --git a/test/test_tree_2d_part2.jl b/test/test_tree_2d_part2.jl index d8e86d14f18..622f12109ff 100644 --- a/test/test_tree_2d_part2.jl +++ b/test/test_tree_2d_part2.jl @@ -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") diff --git a/test/test_unit.jl b/test/test_unit.jl index aa7840ec7aa..20368481607 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -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) From 4fecc09b927bc516258513a673882e6f8b446353 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 17 Nov 2023 09:47:31 +0100 Subject: [PATCH 2/3] set version to v0.6.1 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index da618dc78c1..e324de5fc77 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.6.1-pre" +version = "0.6.1" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 09bdf97df69f55075f853842a97fb2cc2faf703d Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 17 Nov 2023 09:47:49 +0100 Subject: [PATCH 3/3] set development version to v0.6.2-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e324de5fc77..17c02a2adac 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.6.1" +version = "0.6.2-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"