diff --git a/examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl b/examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl
new file mode 100644
index 00000000000..53564b28c4a
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl
@@ -0,0 +1,125 @@
+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
+    return flux(u_inner, normal_direction, equations)
+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)
+
+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
diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl
index 666455fe86c..3f98a3bd59d 100644
--- a/src/equations/compressible_euler_multicomponent_2d.jl
+++ b/src/equations/compressible_euler_multicomponent_2d.jl
@@ -681,6 +681,38 @@ 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::SVector{N,T}, u_rr::SVector{N,T}, normal_direction::SVector{2,T},
+                                     equations::CompressibleEulerMulticomponentEquations2D) where {N, T<:AbstractFloat}
+    # 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 = T(density(u_ll, equations))
+    rho_rr = T(density(u_rr, equations))
+    gamma_ll = T(totalgamma(u_ll, equations))
+    gamma_rr = T(totalgamma(u_rr, equations))
+
+    # Velocity components
+    v_ll_vec = SVector(rho_v1_ll / rho_ll, rho_v2_ll / rho_ll)
+    v_rr_vec = SVector(rho_v1_rr / rho_rr, rho_v2_rr / rho_rr)
+
+    # Project velocities onto the direction normal_direction.
+    v_ll = dot(v_ll_vec, normal_direction)
+    v_rr = dot(v_rr_vec, normal_direction)
+
+    # Compute pressures
+    p_ll = (gamma_ll - one(T)) * (rho_e_ll - T(0.5) * dot(v_ll_vec, v_ll_vec) * rho_ll)
+    p_rr = (gamma_rr - one(T)) * (rho_e_rr - T(0.5) * 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)
+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)
diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl
index 55401a1d129..cbebf6b1a09 100644
--- a/test/test_p4est_2d.jl
+++ b/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,