diff --git a/NEWS.md b/NEWS.md index 4b96e1e2834..0c78484a782 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,6 +12,7 @@ for human readability. - Non-uniform `TreeMesh` available for hyperbolic-parabolic equations. - Capability to set truly discontinuous initial conditions in 1D. - Wetting and drying feature and examples for 1D and 2D shallow water equations +- Implementation of the polytropic Euler equations in 2D - Subcell positivity limiting support for conservative variables in 2D for `TreeMesh` #### Changed diff --git a/examples/structured_2d_dgsem/elixir_eulerpolytropic_convergence.jl b/examples/structured_2d_dgsem/elixir_eulerpolytropic_convergence.jl new file mode 100644 index 00000000000..4fc9281e7a0 --- /dev/null +++ b/examples/structured_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 = flux_hll, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (0.0, 0.0) +coordinates_max = (1.0, 1.0) + +cells_per_dimension = (4, 4) + +mesh = StructuredMesh(cells_per_dimension, + coordinates_min, + coordinates_max) + +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/examples/structured_2d_dgsem/elixir_eulerpolytropic_ec.jl b/examples/structured_2d_dgsem/elixir_eulerpolytropic_ec.jl new file mode 100644 index 00000000000..de15f6b2bc3 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_eulerpolytropic_ec.jl @@ -0,0 +1,80 @@ + +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_weak_blast_wave + +############################################################################### +# Get the DG approximation space + +volume_flux = flux_winters_etal +solver = DGSEM(polydeg=4, surface_flux=flux_winters_etal, + volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the curved quad mesh from a mapping function + +# Mapping as described in https://arxiv.org/abs/2012.12040, but reduced to 2D +function mapping(xi_, eta_) + # Transform input variables between -1 and 1 onto [0,3] + xi = 1.5 * xi_ + 1.5 + eta = 1.5 * eta_ + 1.5 + + y = eta + 3/8 * (cos(1.5 * pi * (2 * xi - 3)/3) * + cos(0.5 * pi * (2 * eta - 3)/3)) + + x = xi + 3/8 * (cos(0.5 * pi * (2 * xi - 3)/3) * + cos(2 * pi * (2 * y - 3)/3)) + + return SVector(x, y) +end + +cells_per_dimension = (16, 16) + +# Create curved mesh with 16 x 16 elements +mesh = StructuredMesh(cells_per_dimension, mapping) + +############################################################################### +# create the semi discretization object + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=100, + save_initial_solution=true, + save_final_solution=true) + +stepsize_callback = StepsizeCallback(cfl=1.0) + +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/examples/structured_2d_dgsem/elixir_eulerpolytropic_isothermal_wave.jl b/examples/structured_2d_dgsem/elixir_eulerpolytropic_isothermal_wave.jl new file mode 100644 index 00000000000..4ab90957579 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_eulerpolytropic_isothermal_wave.jl @@ -0,0 +1,83 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the polytropic Euler equations + +gamma = 1.0 # With gamma = 1 the system is isothermal. +kappa = 1.0 # Scaling factor for the pressure. +equations = PolytropicEulerEquations2D(gamma, kappa) + +# Linear pressure wave in the negative x-direction. +function initial_condition_wave(x, t, equations::PolytropicEulerEquations2D) + rho = 1.0 + v1 = 0.0 + if x[1] > 0.0 + rho = ((1.0 + 0.01 * sin(x[1] * 2 * pi)) / equations.kappa)^(1 / equations.gamma) + v1 = ((0.01 * sin((x[1] - 1 / 2) * 2 * pi)) / equations.kappa) + end + v2 = 0.0 + + return prim2cons(SVector(rho, v1, v2), equations) +end +initial_condition = initial_condition_wave + +volume_flux = flux_winters_etal +solver = DGSEM(polydeg = 2, surface_flux = flux_hll, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-2.0, -1.0) +coordinates_max = (2.0, 1.0) + +cells_per_dimension = (64, 32) + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic) + +mesh = StructuredMesh(cells_per_dimension, + coordinates_min, + coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 200 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 50, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.7) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (1.0e-4, 1.0e-4), + variables = (Trixi.density, pressure)) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, 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); + +# Print the timer summary +summary_callback() diff --git a/examples/structured_2d_dgsem/elixir_eulerpolytropic_wave.jl b/examples/structured_2d_dgsem/elixir_eulerpolytropic_wave.jl new file mode 100644 index 00000000000..fd332b20aef --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_eulerpolytropic_wave.jl @@ -0,0 +1,80 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the polytropic Euler equations + +gamma = 2.0 # Adiabatic monatomic gas in 2d. +kappa = 0.5 # Scaling factor for the pressure. +equations = PolytropicEulerEquations2D(gamma, kappa) + +# Linear pressure wave in the negative x-direction. +function initial_condition_wave(x, t, equations::PolytropicEulerEquations2D) + rho = 1.0 + v1 = 0.0 + if x[1] > 0.0 + rho = ((1.0 + 0.01 * sin(x[1] * 2 * pi)) / equations.kappa)^(1 / equations.gamma) + v1 = ((0.01 * sin((x[1] - 1 / 2) * 2 * pi)) / equations.kappa) + end + v2 = 0.0 + + return prim2cons(SVector(rho, v1, v2), equations) +end +initial_condition = initial_condition_wave + +volume_flux = flux_winters_etal +solver = DGSEM(polydeg = 2, surface_flux = flux_hll, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-2.0, -1.0) +coordinates_max = (2.0, 1.0) + +cells_per_dimension = (64, 32) + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic) + +mesh = StructuredMesh(cells_per_dimension, + coordinates_min, + coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 200 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 50, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.7) + +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); + +# Print the timer summary +summary_callback() diff --git a/src/Trixi.jl b/src/Trixi.jl index b65d03e7975..4169700898a 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -150,7 +150,8 @@ export AcousticPerturbationEquations2D, ShallowWaterEquations1D, ShallowWaterEquations2D, ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D, ShallowWaterEquationsQuasi1D, - LinearizedEulerEquations2D + LinearizedEulerEquations2D, + PolytropicEulerEquations2D export LaplaceDiffusion1D, LaplaceDiffusion2D, CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D, @@ -165,7 +166,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_kennedy_gruber, flux_shima_etal, flux_ec, flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_es_fjordholm_etal, flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal, - flux_chan_etal, flux_nonconservative_chan_etal, + flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal, hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal, # TODO: TrixiShallowWater: move anything with "chen_noelle" to new file hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, diff --git a/src/auxiliary/math.jl b/src/auxiliary/math.jl index 4ecf7dd3fcc..38ea0bda8c8 100644 --- a/src/auxiliary/math.jl +++ b/src/auxiliary/math.jl @@ -211,4 +211,71 @@ Return `x` if `x` is negative, else zero. In other words, return @inline function negative_part(x) return min(x, zero(x)) end + +""" + stolarsky_mean(x, y, gamma) + +Compute an instance of a weighted Stolarsky mean of the form + + stolarsky_mean(x, y, gamma) = (gamma - 1)/gamma * (y^gamma - x^gamma) / (y^(gamma-1) - x^(gamma-1)) + +where `gamma > 1`. + +Problem: The formula above has a removable singularity at `x == y`. Thus, +some care must be taken to implement it correctly without problems or loss +of accuracy when `x ≈ y`. Here, we use the approach proposed by +Winters et al. (2020). +Set f = (y - x) / (y + x) and g = gamma (for compact notation). +Then, we use the expansions + + ((1+f)^g - (1-f)^g) / g = 2*f + (g-1)(g-2)/3 * f^3 + (g-1)(g-2)(g-3)(g-4)/60 * f^5 + O(f^7) + +and + + ((1+f)^(g-1) - (1-f)^(g-1)) / (g-1) = 2*f + (g-2)(g-3)/3 * f^3 + (g-2)(g-3)(g-4)(g-5)/60 * f^5 + O(f^7) + +Inserting the first few terms of these expansions and performing polynomial long division +we find that + + stolarsky_mean(x, y, gamma) ≈ (y + x) / 2 * (1 + (g-2)/3 * f^2 - (g+1)(g-2)(g-3)/45 * f^4 + (g+1)(g-2)(g-3)(2g(g-2)-9)/945 * f^6) + +Since divisions are usually more expensive on modern hardware than +multiplications (Agner Fog), we try to avoid computing two divisions. Thus, +we use + + f^2 = (y - x)^2 / (x + y)^2 + = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y) + +Given ε = 1.0e-4, we use the following algorithm. + + if f^2 < ε + # use the expansion above + else + # use the direct formula (gamma - 1)/gamma * (y^gamma - x^gamma) / (y^(gamma-1) - x^(gamma-1)) + end + +# References +- Andrew R. Winters, Christof Czernik, Moritz B. Schily & Gregor J. Gassner (2020) + Entropy stable numerical approximations for the isothermal and polytropic + Euler equations + [DOI: 10.1007/s10543-019-00789-w](https://doi.org/10.1007/s10543-019-00789-w) +- Agner Fog. + Lists of instruction latencies, throughputs and micro-operation breakdowns + for Intel, AMD, and VIA CPUs. + [https://www.agner.org/optimize/instruction_tables.pdf](https://www.agner.org/optimize/instruction_tables.pdf) +""" +@inline function stolarsky_mean(x, y, gamma) + epsilon_f2 = 1.0e-4 + f2 = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y) # f2 = f^2 + if f2 < epsilon_f2 + # convenience coefficients + c1 = (1 / 3) * (gamma - 2) + c2 = -(1 / 15) * (gamma + 1) * (gamma - 3) * c1 + c3 = -(1 / 21) * (2 * gamma * (gamma - 2) - 9) * c2 + return 0.5 * (x + y) * @evalpoly(f2, 1, c1, c2, c3) + else + return (gamma - 1) / gamma * (y^gamma - x^gamma) / + (y^(gamma - 1) - x^(gamma - 1)) + end +end end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 9bae563d85f..3142dcc2765 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -371,6 +371,11 @@ abstract type AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCO include("compressible_euler_multicomponent_1d.jl") include("compressible_euler_multicomponent_2d.jl") +# PolytropicEulerEquations +abstract type AbstractPolytropicEulerEquations{NDIMS, NVARS} <: + AbstractEquations{NDIMS, NVARS} end +include("polytropic_euler_2d.jl") + # Retrieve number of components from equation instance for the multicomponent case @inline function ncomponents(::AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl new file mode 100644 index 00000000000..d4902bbafb2 --- /dev/null +++ b/src/equations/polytropic_euler_2d.jl @@ -0,0 +1,257 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" + PolytropicEulerEquations2D(gamma, kappa) + +The polytropic Euler equations +```math +\frac{\partial}{\partial t} +\begin{pmatrix} +\rho \\ \rho v_1 \\ \rho v_2 +\end{pmatrix} ++ +\frac{\partial}{\partial x} +\begin{pmatrix} + \rho v_1 \\ \rho v_1^2 + \kappa\rho^\gamma \\ \rho v_1 v_2 +\end{pmatrix} ++ +\frac{\partial}{\partial y} +\begin{pmatrix} +\rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + \kappa\rho^\gamma +\end{pmatrix} += +\begin{pmatrix} +0 \\ 0 \\ 0 +\end{pmatrix} +``` +for an ideal gas with ratio of specific heats `gamma` +in two space dimensions. +Here, ``\rho`` is the density and ``v_1`` and`v_2` the velocities and +```math +p = \kappa\rho^\gamma +``` +the pressure, which we replaced using this relation. +""" +struct PolytropicEulerEquations2D{RealT <: Real} <: + AbstractPolytropicEulerEquations{2, 3} + gamma::RealT # ratio of specific heats + kappa::RealT # fluid scaling factor + + function PolytropicEulerEquations2D(gamma, kappa) + gamma_, kappa_ = promote(gamma, kappa) + new{typeof(gamma_)}(gamma_, kappa_) + end +end + +function varnames(::typeof(cons2cons), ::PolytropicEulerEquations2D) + ("rho", "rho_v1", "rho_v2") +end +varnames(::typeof(cons2prim), ::PolytropicEulerEquations2D) = ("rho", "v1", "v2") + +""" + initial_condition_convergence_test(x, t, equations::PolytropicEulerEquations2D) + +Manufactured smooth initial condition used for convergence tests +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) + + return SVector(h, h / 2, 3 * h / 2) +end + +""" + source_terms_convergence_test(u, x, t, equations::PolytropicEulerEquations2D) + +Source terms used for convergence tests in combination with +[`initial_condition_convergence_test`](@ref). +""" +@inline function source_terms_convergence_test(u, x, t, + equations::PolytropicEulerEquations2D) + 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) + + 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 + + return SVector(r_1, r_2, r_3) +end + +""" + initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquations2D) + +A weak blast wave adapted from +- Sebastian Hennemann, Gregor J. Gassner (2020) + A provably entropy stable subcell shock capturing approach for high order split form DG + [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) +""" +function initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquations2D) + # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3) + # Set up polar coordinates + inicenter = (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) + + # 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) + + return prim2cons(SVector(rho, v1, v2), equations) +end + +# Calculate 1D 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) + rho, v1, v2 = cons2prim(u, equations) + p = pressure(u, equations) + + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + rho_v_normal = rho * v_normal + f1 = rho_v_normal + f2 = rho_v_normal * v1 + p * normal_direction[1] + f3 = rho_v_normal * v2 + p * normal_direction[2] + return SVector(f1, f2, f3) +end + +""" + flux_winters_etal(u_ll, u_rr, normal_direction, + equations::PolytropicEulerEquations2D) + +Entropy conserving two-point flux for isothermal or polytropic gases. +Requires a special weighted Stolarsky mean for the evaluation of the density +denoted here as `stolarsky_mean`. Note, for isothermal gases where `gamma = 1` +this `stolarsky_mean` becomes the [`ln_mean`](@ref). + +For details see Section 3.2 of the following reference +- Andrew R. Winters, Christof Czernik, Moritz B. Schily & Gregor J. Gassner (2020) + Entropy stable numerical approximations for the isothermal and polytropic + Euler equations + [DOI: 10.1007/s10543-019-00789-w](https://doi.org/10.1007/s10543-019-00789-w) +""" +@inline function flux_winters_etal(u_ll, u_rr, normal_direction::AbstractVector, + 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 + 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] + + # 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) + + # Calculate fluxes depending on normal_direction + f1 = rho_mean * 0.5 * (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] + + 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) + 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 + + 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] + + norm_ = norm(normal_direction) + # The v_normals are already scaled by the norm + lambda_min = v_normal_ll - sqrt(equations.gamma * p_ll / rho_ll) * norm_ + lambda_max = v_normal_rr + sqrt(equations.gamma * p_rr / rho_rr) * norm_ + + return lambda_min, lambda_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)) + + return abs(v1) + c, abs(v2) + c +end + +# Convert conservative variables to primitive +@inline function cons2prim(u, equations::PolytropicEulerEquations2D) + rho, rho_v1, rho_v2 = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + + return SVector(rho, v1, v2) +end + +# Convert conservative variables to entropy +@inline function cons2entropy(u, equations::PolytropicEulerEquations2D) + rho, rho_v1, rho_v2 = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + 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 + internal_energy = equations.kappa * log(rho) + else # equations.gamma > 1 # polytropic gas + internal_energy = equations.kappa * rho^(equations.gamma - 1) / + (equations.gamma - 1.0) + end + + w1 = internal_energy + p / rho - 0.5 * v_square + w2 = v1 + w3 = v2 + + return SVector(w1, w2, w3) +end + +# Convert primitive to conservative variables +@inline function prim2cons(prim, equations::PolytropicEulerEquations2D) + rho, v1, v2 = prim + rho_v1 = rho * v1 + rho_v2 = rho * v2 + return SVector(rho, rho_v1, rho_v2) +end + +@inline function density(u, equations::PolytropicEulerEquations2D) + rho = u[1] + return rho +end + +@inline function pressure(u, equations::PolytropicEulerEquations2D) + rho, rho_v1, rho_v2 = u + p = equations.kappa * rho^equations.gamma + return p +end +end # @muladd diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 75937ba82ad..6d528abc7af 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -220,6 +220,30 @@ isdir(outdir) && rm(outdir, recursive=true) tspan = (0.0, 0.3)) end + @trixi_testset "elixir_eulerpolytropic_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"), + l2 = [0.0016688820596537988, 0.0025921681885685425, 0.003280950351435014], + linf = [0.010994679664394269, 0.01331197845637, 0.020080117011346488]) + end + + @trixi_testset "elixir_eulerpolytropic_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_ec.jl"), + l2 = [0.03647890611450939, 0.025284915444045052, 0.025340697771609126], + linf = [0.32516731565355583, 0.37509762516540046, 0.29812843284727336]) + end + + @trixi_testset "elixir_eulerpolytropic_isothermal_wave.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_isothermal_wave.jl"), + l2 = [0.004998778491726366, 0.004998916000294425, 9.259136963058664e-17], + linf = [0.010001103673834888, 0.010051165098399503, 7.623942913643681e-16]) + end + + @trixi_testset "elixir_eulerpolytropic_wave.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_wave.jl"), + l2 = [0.23642682112204072, 0.20904264390331334, 8.174982691297391e-17], + linf = [0.4848250368349989, 0.253350873815695, 4.984552457753618e-16]) + end + @trixi_testset "elixir_hypdiff_nonperiodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_nonperiodic.jl"), l2 = [0.8799744480157664, 0.8535008397034816, 0.7851383019164209],