Added abstract vector for max_abs_speed_naive and CompressibleEulerMulticomponentEquations2D.

examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl

@@ -0,0 +1,127 @@
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+
+# 1) Dry Air  2) SF_6
+equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648),
+                                                       gas_constants = (0.287, 1.578))
+
+"""
+    initial_condition_shock(coordinates, t, equations::CompressibleEulerEquations2D)
+
+Shock traveling from left to right where it interacts with a Perturbed interface.
+"""
+@inline function initial_condition_shock(x, t,
+                                        equations::CompressibleEulerMulticomponentEquations2D)
+    rho_0 = 1.25 # kg/m^3
+    p_0 = 101325 # Pa
+    T_0 = 293 # K
+    u_0 = 352 # m/s
+    d = 20
+    w = 40
+
+    if x[1] < 25
+        # Shock region.
+        v1 = 0.35
+        v2 = 0.0
+        rho1 = 1.72
+        rho2 = 0.03
+        p = 1.57
+    elseif (x[1] <= 30) || (x[1] <= 70 && abs(125 - x[2]) > w/2)
+        # Intermediate region.
+        v1 = 0.0
+        v2 = 0.0
+        rho1 = 1.25
+        rho2 = 0.03
+        p = 1.0
+    else (x[1] <= 70 + d)
+        # SF_6 region.
+        v1 = 0.0
+        v2 = 0.0
+        rho1 = 0.03
+        rho2 = 6.03 #SF_6
+        p = 1.0
+    end
+
+    return prim2cons(SVector(v1, v2, p, rho1, rho2), equations)
+end
+
+# Define the simulation.
+initial_condition = initial_condition_shock
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+basis = LobattoLegendreBasis(3)
+
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                positivity_variables_cons = ["rho" * string(i)
+                                                             for i in eachcomponent(equations)])
+
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (250.0, 250.0)
+mesh = P4estMesh((32, 32), polydeg=3, coordinates_min=(0.0, 0.0),
+                 coordinates_max=(250.0, 250.0), periodicity=(false, true),
+                 initial_refinement_level=0)
+
+# Completely free outflow
+function boundary_condition_outflow(u_inner, normal_direction::AbstractVector,
+                                    x, t,
+                                    surface_flux_function,
+                                    equations)
+    # Calculate the boundary flux entirely from the internal solution state
+    flux = Trixi.flux(u_inner, normal_direction, equations)
+
+    return flux
+end
+
+boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition),
+                           :x_pos => boundary_condition_outflow)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 10
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = true)
+#                                      extra_analysis_integrals = (enstrophy_multi_euler,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(dt = 2.0,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.2)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+stage_callbacks = (SubcellLimiterIDPCorrection(),
+                   BoundsCheckCallback(save_errors = false, interval = 100))
+# `interval` is used when calling this elixir in the tests with `save_errors=true`.
+
+@time begin
+    sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
+                    dt = 0.1, # solve needs some value here but it will be overwritten by the stepsize_callback
+                    ode_default_options()..., callback = callbacks);
+end

src/equations/compressible_euler_multicomponent_2d.jl

@@ -681,6 +681,40 @@ end
     return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr)
 end
 
+# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
+@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::CompressibleEulerMulticomponentEquations2D)
+    # Unpack conservative variables
+    rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
+    rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr
+
+    # Get densities and gammas
+    rho_ll = density(u_ll, equations)
+    rho_rr = density(u_rr, equations)
+    gamma_ll = totalgamma(u_ll, equations)
+    gamma_rr = totalgamma(u_rr, equations)
+
+    # Velocity components
+    v_ll_vec = SVector(rho_v1_ll, rho_v2_ll) ./ rho_ll
+    v_rr_vec = SVector(rho_v1_rr, rho_v2_rr) ./ rho_rr
+    # Normalize the direction vector
+    norm_ = norm(normal_direction)
+    normal_vector = normal_direction / norm_
+    # Project velocities onto the direction normal_direction.
+    v_ll = dot(v_ll_vec, normal_vector)
+    v_rr = dot(v_rr_vec, normal_vector)
+
+    # Compute pressures
+    p_ll = (gamma_ll - 1) * (rho_e_ll - 0.5f0 * dot(v_ll_vec, v_ll_vec) * rho_ll)
+    p_rr = (gamma_rr - 1) * (rho_e_rr - 0.5f0 * dot(v_rr_vec, v_rr_vec) * rho_rr)
+
+    # Sound speeds
+    c_ll = sqrt(gamma_ll * p_ll / rho_ll)
+    c_rr = sqrt(gamma_rr * p_rr / rho_rr)
+
+    return (max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr)) * norm_
+end
+
 # Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
 @inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
                                equations::CompressibleEulerMulticomponentEquations2D)

test/test_p4est_2d.jl

@@ -560,6 +560,27 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "elixir_eulermulti_shock.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock.jl"),
+                        l2=[
+                            0.09160388140703994,
+                            0.000736779618610053,
+                            0.22564799038297123,
+                            0.09074993408575387,
+                            0.2714279606670099
+                        ],
+                        linf=[
+                            0.7071145839509081,
+                            0.031436413121202246,
+                            1.7750084462888234,
+                            1.1624061471508902,
+                            4.335595275270981
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
 @trixi_testset "elixir_mhd_alfven_wave.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
                         l2=[1.0513414461545583e-5, 1.0517900957166411e-6,