@testsnippet P4estMesh3D begin
    EXAMPLES_DIR = joinpath(examples_dir(), "p4est_3d_dgsem")
end

@testitem "P4estMesh3D: elixir_advection_basic.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
                        # Expected errors are exactly the same as with TreeMesh!
                        l2=[0.00016263963870641478],
                        linf=[0.0014537194925779984])
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_advection_limiter_liu_zhang.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem",
                                 "elixir_advection_limiter_liu_zhang.jl"),
                        mesh=P4estMesh((4, 4, 4), polydeg = 3,
                                       coordinates_min = (-1.0, -1.0, -1.0),
                                       coordinates_max = (1.0, 1.0, 1.0),
                                       initial_refinement_level = 1,
                                       periodicity = true),
                        l2=[0.3917747509079011],
                        linf=[1.3569114615249707])
    u = Trixi.wrap_array_native(sol.u[end], semi)
    # matches thresholds = (1e-3,) up to a tolerance
    @test minimum(u) > 1e-3 - 10 * eps()
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_advection_unstructured_curved.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_advection_unstructured_curved.jl"),
                        l2=[0.0004750004258546538],
                        linf=[0.026527551737137167])
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_advection_nonconforming.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming.jl"),
                        l2=[0.00253595715323843],
                        linf=[0.016486952252155795])
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_advection_amr.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr.jl"),
                        # Expected errors are exactly the same as with TreeMesh!
                        l2=[9.773852895157622e-6],
                        linf=[0.0005853874124926162])
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_advection_amr_unstructured_curved.jl" setup=[
    Setup,
    P4estMesh3D
] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_advection_amr_unstructured_curved.jl"),
                        l2=[1.6163120948209677e-5],
                        linf=[0.0010572201890564834],
                        tspan=(0.0, 1.0))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_advection_cubed_sphere.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_cubed_sphere.jl"),
                        l2=[0.002006918015656413],
                        linf=[0.027655117058380085])
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_advection_restart.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
                        l2=[0.002590388934758452],
                        linf=[0.01840757696885409])
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_source_terms_nonconforming_unstructured_curved.jl" setup=[
    Setup,
    P4estMesh3D
] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_euler_source_terms_nonconforming_unstructured_curved.jl"),
                        l2=[
                            4.070355207909268e-5,
                            4.4993257426833716e-5,
                            5.10588457841744e-5,
                            5.102840924036687e-5,
                            0.00019986264001630542
                        ],
                        linf=[
                            0.0016987332417202072,
                            0.003622956808262634,
                            0.002029576258317789,
                            0.0024206977281964193,
                            0.008526972236273522
                        ],
                        tspan=(0.0, 0.01))

    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_euler_source_terms_nonconforming_unstructured_curved.jl"),
                        surface_flux=FluxPlusDissipation(flux_ranocha,
                                                         DissipationMatrixWintersEtal()),
                        l2=[
                            4.068002997087932e-5,
                            4.4742882348806466e-5,
                            5.101817697733163e-5,
                            5.100410876233901e-5,
                            0.000199848133462063
                        ],
                        linf=[
                            0.0013080357114820806,
                            0.0028524316301083985,
                            0.0019100643150573582,
                            0.0024800222220195955,
                            0.00830424488849335
                        ],
                        tspan=(0.0, 0.01))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_source_terms_nonperiodic.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_euler_source_terms_nonperiodic.jl"),
                        l2=[
                            0.0015106060984283647,
                            0.0014733349038567685,
                            0.00147333490385685,
                            0.001473334903856929,
                            0.0028149479453087093
                        ],
                        linf=[
                            0.008070806335238156,
                            0.009007245083113125,
                            0.009007245083121784,
                            0.009007245083102688,
                            0.01562861968368434
                        ],
                        tspan=(0.0, 1.0))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_free_stream.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
                        l2=[
                            5.162664597942288e-15,
                            1.941857343642486e-14,
                            2.0232366394187278e-14,
                            2.3381518645408552e-14,
                            7.083114561232324e-14
                        ],
                        linf=[
                            7.269740365245525e-13,
                            3.289868377720495e-12,
                            4.440087186807773e-12,
                            3.8686831516088205e-12,
                            9.412914891981927e-12
                        ],
                        tspan=(0.0, 0.03))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_free_stream.jl (VolumeIntegralPureLGLFiniteVolume)" setup=[
    Setup,
    P4estMesh3D
] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
                        volume_integral=VolumeIntegralPureLGLFiniteVolume(flux_hllc),
                        cfl=1.0,
                        l2=[
                            1.9974155101075348e-15,
                            1.1961083306610996e-14,
                            1.2872974022570237e-14,
                            1.5462968285712546e-14,
                            2.6094715468434938e-14
                        ],
                        linf=[
                            2.828848266744899e-13,
                            1.6810441927361808e-12,
                            2.784633634789202e-12,
                            2.8425040099477883e-12,
                            3.7339020764193265e-12
                        ],
                        tspan=(0.0, 0.03))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_free_stream_extruded.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_extruded.jl"),
                        l2=[
                            8.444868392439035e-16,
                            4.889826056731442e-15,
                            2.2921260987087585e-15,
                            4.268460455702414e-15,
                            1.1356712092620279e-14
                        ],
                        linf=[
                            7.749356711883593e-14,
                            2.8792246364872653e-13,
                            1.1121659149182506e-13,
                            3.3228975127030935e-13,
                            9.592326932761353e-13
                        ],
                        tspan=(0.0, 0.1))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_free_stream_extruded_fvO2.jl" setup=[
    Setup,
    P4estMesh3D
] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_euler_free_stream_extruded_fvO2.jl"),
                        l2=[
                            1.5286679406524262e-16,
                            4.950311927136664e-16,
                            4.69913693630953e-16,
                            4.012046698403924e-16,
                            2.0419659904946413e-15
                        ],
                        linf=[
                            2.220446049250313e-15,
                            3.3029134982598407e-15,
                            7.716050021144838e-15,
                            1.27675647831893e-14,
                            3.019806626980426e-14
                        ],
                        tspan=(0.0, 0.1))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_free_stream_boundaries.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_euler_free_stream_boundaries.jl"),
                        l2=[
                            6.530157034651212e-16, 1.6057829680004379e-15,
                            3.31107455378537e-15, 3.908829498281281e-15,
                            5.048390610424672e-15
                        ],
                        linf=[
                            4.884981308350689e-15, 1.1921019726912618e-14,
                            1.5432100042289676e-14, 2.298161660974074e-14,
                            6.039613253960852e-14
                        ])
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_free_stream_boundaries_float32.jl" setup=[
    Setup,
    P4estMesh3D
] tags=[:p4est_part2] begin
    # Expected errors are taken from elixir_euler_free_stream_boundaries.jl
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_euler_free_stream_boundaries_float32.jl"),
                        l2=[
                            Float32(6.530157034651212e-16),
                            Float32(1.6057829680004379e-15),
                            Float32(3.31107455378537e-15),
                            Float32(3.908829498281281e-15),
                            Float32(5.048390610424672e-15)
                        ],
                        linf=[
                            Float32(4.884981308350689e-15),
                            Float32(1.1921019726912618e-14),
                            Float32(1.5432100042289676e-14),
                            Float32(2.298161660974074e-14),
                            Float32(6.039613253960852e-14)
                        ],
                        RealT_for_test_tolerances=Float32)
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_free_stream_extruded.jl with HLLC FLux" setup=[
    Setup,
    P4estMesh3D
] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_extruded.jl"),
                        l2=[
                            8.444868392439035e-16,
                            4.889826056731442e-15,
                            2.2921260987087585e-15,
                            4.268460455702414e-15,
                            1.1356712092620279e-14
                        ],
                        linf=[
                            7.749356711883593e-14,
                            4.513472928735496e-13,
                            2.9790059308254513e-13,
                            1.057154364048074e-12,
                            1.6271428648906294e-12
                        ],
                        tspan=(0.0, 0.1),
                        surface_flux=flux_hllc)
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_ec.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
                        l2=[
                            0.010380390326164493,
                            0.006192950051354618,
                            0.005970674274073704,
                            0.005965831290564327,
                            0.02628875593094754
                        ],
                        linf=[
                            0.3326911600075694,
                            0.2824952141320467,
                            0.41401037398065543,
                            0.45574161423218573,
                            0.8099577682187109
                        ],
                        tspan=(0.0, 0.2))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_ec.jl (flux_chandrashekar)" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
                        l2=[
                            0.010368548525287055,
                            0.006216054794583285,
                            0.006020401857347216,
                            0.006019175682769779,
                            0.026228080232814154
                        ],
                        linf=[
                            0.3169376449662026,
                            0.28950510175646726,
                            0.4402523227566396,
                            0.4869168122387365,
                            0.7999141641954051
                        ],
                        tspan=(0.0, 0.2),
                        volume_flux=flux_chandrashekar)
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_sedov.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"),
                        l2=[
                            7.82070951e-02,
                            4.33260474e-02,
                            4.33260474e-02,
                            4.33260474e-02,
                            3.75260911e-01
                        ],
                        linf=[
                            7.45329845e-01,
                            3.21754792e-01,
                            3.21754792e-01,
                            3.21754792e-01,
                            4.76151527e+00
                        ],
                        tspan=(0.0, 0.3))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_sedov_sc_subcell.jl (positivity bounds)" setup=[
    Setup,
    P4estMesh3D
] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_sc_subcell.jl"),
                        l2=[
                            0.19427117014566905,
                            0.07557679851688888,
                            0.07557679851688896,
                            0.07557679851688917,
                            0.3713893373335463
                        ],
                        linf=[
                            2.754292081408987,
                            1.888627709320322,
                            1.8886277093203232,
                            1.888627709320322,
                            4.971280431903264
                        ],
                        tspan=(0.0, 0.3))
    limiter = semi.solver.volume_integral.limiter
    deviations = collect(values(limiter.cache.idp_bounds_delta_global))
    @test all(isfinite, deviations)
    @test maximum(deviations) <= 1.0e-13

    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    # Larger values for allowed allocations due to usage of custom
    # integrator which are not *recorded* for the methods from
    # OrdinaryDiffEq.jl
    # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 15_000)
end

@testitem "P4estMesh3D: elixir_euler_sedov_sc_subcell.jl (local bounds)" setup=[
    Setup,
    P4estMesh3D
] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_sc_subcell.jl"),
                        local_twosided_variables_cons=["rho"],
                        local_onesided_variables_nonlinear=[(entropy_guermond_etal,
                                                             min)],
                        max_iterations_newton=30,
                        l2=[
                            0.16504564013491585,
                            0.06461384162458203,
                            0.06461384162461223,
                            0.06461384162461678,
                            0.36193245790622036
                        ],
                        linf=[
                            0.9138327077620716,
                            0.5707102472596818,
                            0.5707102472739252,
                            0.5707102472781822,
                            4.777595503303726
                        ],
                        tspan=(0.0, 0.3))
    limiter = semi.solver.volume_integral.limiter
    deviations = collect(values(limiter.cache.idp_bounds_delta_global))
    @test all(isfinite, deviations)
    @test maximum(deviations) <= 1.0e-13

    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    # Larger values for allowed allocations due to usage of custom
    # integrator which are not *recorded* for the methods from
    # OrdinaryDiffEq.jl
    # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 15_000)
end

@testitem "P4estMesh3D: elixir_euler_sedov_sc_subcell.jl (local bounds, nonperiodic)" setup=[
    Setup,
    P4estMesh3D
] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_sc_subcell.jl"),
                        local_twosided_variables_cons=["rho"],
                        local_onesided_variables_nonlinear=[(entropy_guermond_etal,
                                                             min)],
                        max_iterations_newton=50,
                        periodicity=false,
                        boundary_conditions=BoundaryConditionDirichlet(initial_condition),
                        l2=[
                            0.16504564013491585,
                            0.06461384162458203,
                            0.06461384162461223,
                            0.06461384162461678,
                            0.36193245790622036
                        ],
                        linf=[
                            0.9138327077620716,
                            0.5707102472596818,
                            0.5707102472739252,
                            0.5707102472781822,
                            4.777595503303726
                        ],
                        tspan=(0.0, 0.3))
    limiter = semi.solver.volume_integral.limiter
    deviations = collect(values(limiter.cache.idp_bounds_delta_global))
    @test all(isfinite, deviations)
    @test maximum(deviations) <= 1.0e-13

    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    # Larger values for allowed allocations due to usage of custom
    # integrator which are not *recorded* for the methods from
    # OrdinaryDiffEq.jl
    # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 15_000)
end

@testitem "P4estMesh3D: elixir_euler_sedov.jl (HLLE)" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"),
                        l2=[
                            0.09946224487902565,
                            0.04863386374672001,
                            0.048633863746720116,
                            0.04863386374672032,
                            0.3751015774232693
                        ],
                        linf=[
                            0.789241521871487,
                            0.42046970270100276,
                            0.42046970270100276,
                            0.4204697027010028,
                            4.730877375538398
                        ],
                        tspan=(0.0, 0.3),
                        surface_flux=flux_hlle)
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_source_terms_nonconforming_earth.jl" setup=[
    Setup,
    P4estMesh3D
] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_euler_source_terms_nonconforming_earth.jl"),
                        l2=[
                            6.040180337738628e-6,
                            5.4254175153621895e-6,
                            5.677698851333843e-6,
                            5.8017136892469794e-6,
                            1.3637854615117974e-5
                        ],
                        linf=[
                            0.00013996924184311865,
                            0.00013681539559939893,
                            0.00013681539539733834,
                            0.00013681539541021692,
                            0.00016833038543762058
                        ],
                        # Decrease tolerance of adaptive time stepping to get similar results across different systems
                        abstol=1.0e-11, reltol=1.0e-11)
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_circular_wind_nonconforming.jl" setup=[
    Setup,
    P4estMesh3D
] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_euler_circular_wind_nonconforming.jl"),
                        l2=[
                            1.5737711609657832e-7,
                            3.8630261900166194e-5,
                            3.8672287531936816e-5,
                            3.6865116098660796e-5,
                            0.05508620970403884
                        ],
                        linf=[
                            2.268845333053271e-6,
                            0.000531462302113539,
                            0.0005314624461298934,
                            0.0005129931254772464,
                            0.7942778058932163
                        ],
                        tspan=(0.0, 2e2))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_baroclinic_instability.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_euler_baroclinic_instability.jl"),
                        l2=[
                            6.725093801700048e-7,
                            0.00021710076010951073,
                            0.0004386796338203878,
                            0.00020836270267103122,
                            0.07601887903440395
                        ],
                        linf=[
                            1.9107530539574924e-5,
                            0.02980358831035801,
                            0.048476331898047564,
                            0.02200137344113612,
                            4.848310144356219
                        ],
                        tspan=(0.0, 1e2),
                        # Decrease tolerance of adaptive time stepping to get similar results across different systems
                        abstol=1.0e-9, reltol=1.0e-9)
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_source_terms_nonperiodic_hohqmesh.jl" setup=[
    Setup,
    P4estMesh3D
] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_euler_source_terms_nonperiodic_hohqmesh.jl"),
                        l2=[
                            0.0042023406458005464,
                            0.004122532789279737,
                            0.0042448149597303616,
                            0.0036361316700401765,
                            0.007389845952982495
                        ],
                        linf=[
                            0.04530610539892499,
                            0.02765695110527666,
                            0.05670295599308606,
                            0.048396544302230504,
                            0.1154589758186293
                        ])
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_tandem_spheres.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_euler_tandem_spheres.jl"),
                        # Strict tols to avoid issues with different systems (CI vs local machines)
                        tspan=(0.0, 1e-3), tols=1e-8,
                        l2=[
                            1.4563595831943405e-6,
                            1.1433846164578946e-6,
                            6.45964041497078e-7,
                            6.446228703680879e-7,
                            3.648646693902689e-6
                        ],
                        linf=[
                            0.08928818652777637,
                            0.09007801634411498,
                            0.04241320785718675,
                            0.04253018852293939,
                            0.22552627764127076
                        ])
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_free_stream_hybrid_mesh.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_euler_free_stream_hybrid_mesh.jl"),
                        l2=[
                            6.459343647140878e-16,
                            1.6402804603545678e-15,
                            2.0956682982886904e-15,
                            2.3267454534698676e-15,
                            3.802615222445529e-15
                        ],
                        linf=[
                            1.1657341758564144e-14,
                            1.6056600493641326e-14,
                            3.5665914666083154e-14,
                            3.68594044175552e-14,
                            6.217248937900877e-14
                        ])
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_mhd_alfven_wave_er.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_mhd_alfven_wave_er.jl"),
                        l2=[
                            0.0052864046546744065,
                            0.009963357787771665,
                            0.006635699953141596,
                            0.01295540589311982,
                            0.013939326496053958,
                            0.010192741315114568,
                            0.004631666336074305,
                            0.012267586777052244,
                            0.0018063823439272181
                        ],
                        linf=[
                            0.021741826900806394,
                            0.0470226920658848,
                            0.025036937229995254,
                            0.05043002191230382,
                            0.06018360063552164,
                            0.04338351710391075,
                            0.023607975939848536,
                            0.050740527490335,
                            0.006909064342577296
                        ])
    # Larger values for allowed allocations due to usage of custom
    # integrator which are not *recorded* for the methods from
    # OrdinaryDiffEq.jl
    # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 15_000)
end

@testitem "P4estMesh3D: elixir_mhd_alfven_wave_nonconforming.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_mhd_alfven_wave_nonconforming.jl"),
                        l2=[
                            0.0001788543743594658,
                            0.000624334205581902,
                            0.00022892869974368887,
                            0.0007223464581156573,
                            0.0006651366626523314,
                            0.0006287275014743352,
                            0.000344484339916008,
                            0.0007179788287557142,
                            8.632896980651243e-7
                        ],
                        linf=[
                            0.0010730565632763867,
                            0.004596749809344033,
                            0.0013235269262853733,
                            0.00468874234888117,
                            0.004719267084104306,
                            0.004228339352211896,
                            0.0037503625505571625,
                            0.005104176909383168,
                            9.738081186490818e-6
                        ],
                        tspan=(0.0, 0.25))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_mhd_alfven_wave_nonperiodic.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_mhd_alfven_wave_nonperiodic.jl"),
                        l2=[
                            0.00017912812934894293,
                            0.000630910737693146,
                            0.0002256138768371346,
                            0.0007301686017397987,
                            0.0006647296256552257,
                            0.0006409790941359089,
                            0.00033986873316986315,
                            0.0007277161123570452,
                            1.3184121257198033e-5
                        ],
                        linf=[
                            0.0012248374096375247,
                            0.004857541490859554,
                            0.001813452620706816,
                            0.004803571938364726,
                            0.005271403957646026,
                            0.004571200760744465,
                            0.002618188297242474,
                            0.005010126350015381,
                            6.309149507784953e-5
                        ],
                        tspan=(0.0, 0.25))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_mhd_shockcapturing_amr.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_amr.jl"),
                        l2=[
                            0.0062973565893792004,
                            0.006436273914579104,
                            0.007112703307027178,
                            0.006529650167358523,
                            0.020607452343745017,
                            0.005560993001492338,
                            0.007576418168749763,
                            0.0055721349394598635,
                            3.8269125984310296e-6
                        ],
                        linf=[
                            0.2090718196650192,
                            0.1863884052971854,
                            0.23475479927204168,
                            0.19460789763442982,
                            0.6859816363887359,
                            0.15171474186273914,
                            0.22404690260234983,
                            0.16808957604979002,
                            0.0005083795485317637
                        ],
                        tspan=(0.0, 0.04))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_mhd_shockcapturing_subcell.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_subcell.jl"),
                        l2=[
                            0.006729931970167595,
                            0.008638393158639436,
                            0.008978257148101689,
                            0.0085466685190268,
                            0.0285664608641833,
                            0.005796806835751598,
                            0.007485378539184046,
                            0.005846216235895686,
                            1.116115482859488e-5
                        ],
                        linf=[
                            0.31507940417652525,
                            0.27581230560179737,
                            0.5096527712168957,
                            0.2900021150087706,
                            0.9484970527977867,
                            0.2591599747174065,
                            0.22934145164154485,
                            0.2868673643088755,
                            0.0013663401454538622
                        ],
                        tspan=(0.0, 0.04))
    limiter = semi.solver.volume_integral.limiter
    deviations = collect(values(limiter.cache.idp_bounds_delta_global))
    @test all(isfinite, deviations)
    @test maximum(deviations) <= 1.0e-13

    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    # Larger values for allowed allocations due to usage of custom
    # integrator which are not *recorded* for the methods from
    # OrdinaryDiffEq.jl
    # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 15_000)
end

@testitem "P4estMesh3D: elixir_mhd_amr_entropy_bounded.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_amr_entropy_bounded.jl"),
                        l2=[
                            0.005430006338127661,
                            0.006186402899876596,
                            0.012171513410597289,
                            0.006181479343504159,
                            0.035068817354117605,
                            0.004967715666538709,
                            0.006592173316509503,
                            0.0050151140388451105,
                            5.146547644807638e-6
                        ],
                        linf=[
                            0.18655204102670386,
                            0.20397573777286138,
                            0.3700839435299759,
                            0.23329319876321034,
                            1.0348619438460904,
                            0.18462694496595722,
                            0.20648634653698617,
                            0.18947822281424997,
                            0.0005083794158781671
                        ],
                        tspan=(0.0, 0.04))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_linearizedeuler_convergence.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_linearizedeuler_convergence.jl"),
                        l2=[
                            0.04452389418193219, 0.03688186699434862,
                            0.03688186699434861, 0.03688186699434858,
                            0.044523894181932186
                        ],
                        linf=[
                            0.2295447498696467, 0.058369658071546704,
                            0.05836965807154648, 0.05836965807154648, 0.2295447498696468
                        ])
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

@testitem "P4estMesh3D: elixir_euler_weak_blast_wave_amr.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"),
                        l2=[
                            0.01374649869395016,
                            0.01993458602992416,
                            0.020403214655756098,
                            0.020408263790073853,
                            0.14975849029503904
                        ],
                        linf=[
                            0.4411601724293266,
                            0.668308654218055,
                            0.7351134068050753,
                            0.5955002383710662,
                            3.1811162616598985
                        ],
                        tspan=(0.0, 0.025))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
    # Check for conservation
    state_integrals = Trixi.integrate(sol.u[2], semi)
    initial_state_integrals = analysis_callback.affect!.initial_state_integrals

    @test isapprox(state_integrals[1], initial_state_integrals[1], atol = 1e-13)
    @test isapprox(state_integrals[2], initial_state_integrals[2], atol = 1e-13)
    @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13)
    @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13)
    @test isapprox(state_integrals[5], initial_state_integrals[5], atol = 1e-13)
end

@testitem "P4estMesh3D: elixir_euler_ONERA_M6_wing.jl" setup=[Setup, P4estMesh3D] tags=[:p4est_part2] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR,
                                 "elixir_euler_ONERA_M6_wing.jl"),
                        l2=[
                            1.3302852203314697e-7,
                            7.016342225152883e-8,
                            1.0954098970860626e-7,
                            6.834890433113107e-8,
                            3.796737956937651e-7
                        ],
                        linf=[
                            0.08856648749331164,
                            0.07431651477033197,
                            0.08791247483932041,
                            0.012973811024139751,
                            0.25575828277482016
                        ],
                        tspan=(0.0, 5e-8))
    # Ensure that we do not have excessive memory allocations
    # (e.g., from type instabilities)
    @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000)
end

# Multi-ion MHD tests
