@testsnippet UnitTests begin
    using LinearAlgebra: norm, dot
    using SparseArrays
    using DelimitedFiles: readdlm
    using ForwardDiff
    using Convex: Convex
    using ECOS: Optimizer
    using NLsolve: nlsolve
    import SparseConnectivityTracer: TracerSparsityDetector, jacobian_eltype,
                                     jacobian_sparsity
    import SparseMatrixColorings: ColoringProblem, GreedyColoringAlgorithm, coloring,
                                  column_colors
    import FiniteDiff: finite_difference_jacobian!
end

@testitem "Unit: Spectral analysis" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    @testset "compute_kinetic_energy_spectrum" begin
        rho_2d = ones(4, 4)
        velocity_1_2d = ones(4, 4)
        velocity_2_2d = zeros(4, 4)

        wavenumbers_2d, energy_spectrum_2d = Trixi.compute_kinetic_energy_spectrum(velocity_1_2d,
                                                                                   velocity_2_2d)
        @test wavenumbers_2d == 0:3
        @test energy_spectrum_2d[1] ≈ 0.5
        @test all(isapprox.(energy_spectrum_2d[2:end], 0, atol = 100 * eps()))

        velocity_1_2d .= [sin(2 * pi * (i - 1) / size(rho_2d, 1))
                          for i in axes(rho_2d, 1), j in axes(rho_2d, 2)]
        velocity_2_2d .= [cos(2 * pi * (j - 1) / size(rho_2d, 2))
                          for i in axes(rho_2d, 1), j in axes(rho_2d, 2)]
        _, energy_spectrum_2d = Trixi.compute_kinetic_energy_spectrum(sqrt.(rho_2d) .*
                                                                      velocity_1_2d,
                                                                      sqrt.(rho_2d) .*
                                                                      velocity_2_2d)
        mean_kinetic_energy_2d = sum(@. 0.5 * rho_2d *
                                        (velocity_1_2d^2 + velocity_2_2d^2)) /
                                 length(rho_2d)
        @test sum(energy_spectrum_2d) ≈ mean_kinetic_energy_2d

        rho_3d = ones(4, 4, 4)
        velocity_1_3d = ones(4, 4, 4)
        velocity_2_3d = zeros(4, 4, 4)
        velocity_3_3d = zeros(4, 4, 4)

        wavenumbers_3d, energy_spectrum_3d = Trixi.compute_kinetic_energy_spectrum(velocity_1_3d,
                                                                                   velocity_2_3d,
                                                                                   velocity_3_3d)
        @test wavenumbers_3d == 0:3
        @test energy_spectrum_3d[1] ≈ 0.5
        @test all(isapprox.(energy_spectrum_3d[2:end], 0, atol = 100 * eps()))

        velocity_1_3d .= [sin(2 * pi * (i - 1) / size(rho_3d, 1))
                          for i in axes(rho_3d, 1), j in axes(rho_3d, 2),
                              k in axes(rho_3d, 3)]
        velocity_2_3d .= [cos(2 * pi * (j - 1) / size(rho_3d, 2))
                          for i in axes(rho_3d, 1), j in axes(rho_3d, 2),
                              k in axes(rho_3d, 3)]
        velocity_3_3d .= [sin(2 * pi * (k - 1) / size(rho_3d, 3))
                          for i in axes(rho_3d, 1), j in axes(rho_3d, 2),
                              k in axes(rho_3d, 3)]
        _, energy_spectrum_3d = Trixi.compute_kinetic_energy_spectrum(sqrt.(rho_3d) .*
                                                                      velocity_1_3d,
                                                                      sqrt.(rho_3d) .*
                                                                      velocity_2_3d,
                                                                      sqrt.(rho_3d) .*
                                                                      velocity_3_3d)
        mean_kinetic_energy_3d = sum(@. 0.5 * rho_3d *
                                        (velocity_1_3d^2 + velocity_2_3d^2 +
                                         velocity_3_3d^2)) /
                                 length(rho_3d)
        @test sum(energy_spectrum_3d) ≈ mean_kinetic_energy_3d
    end
end

@testitem "Unit: SerialTree" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    @testset "constructors" begin
        @test_nowarn Trixi.SerialTree(Val(1), 10, 0.0, 1.0, true)
        @test_nowarn Trixi.SerialTree{1}(10, 0.0, 1.0, true)
    end

    @testset "helper functions" begin
        t = Trixi.SerialTree(Val(1), 10, 0.0, 1.0, true)
        @test_nowarn display(t)
        @test Trixi.ndims(t) == 1
        @test Trixi.has_any_neighbor(t, 1, 1) == true
        @test Trixi.isperiodic(t, 1) == true
        @test Trixi.n_children_per_cell(t) == 2
        @test Trixi.n_directions(t) == 2
    end

    @testset "refine!/coarsen!" begin
        t = Trixi.SerialTree(Val(1), 10, 0.0, 1.0, true)
        @test Trixi.refine!(t) == [1]
        @test Trixi.coarsen!(t) == [1]
        @test Trixi.refine!(t) == [1]
        @test Trixi.coarsen!(t, 1) == [1]
        @test Trixi.coarsen!(t) == Int[] # Coarsen twice to check degenerate case of single-cell tree
        @test Trixi.refine!(t) == [1]
        @test Trixi.refine!(t) == [2, 3]
        @test Trixi.coarsen_box!(t, [-0.5], [0.0]) == [2]
        @test Trixi.coarsen_box!(t, 0.0, 0.5) == [3]
        @test isnothing(Trixi.reset_data_structures!(t))
    end
end

@testitem "Unit: ParallelTree" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    @testset "constructors" begin
        @test_nowarn Trixi.ParallelTree(Val(1), 10, 0.0, 1.0, true)
        @test_nowarn Trixi.ParallelTree{1}(10, 0.0, 1.0, true)
    end

    @testset "helper functions" begin
        t = Trixi.ParallelTree(Val(1), 10, 0.0, 1.0, true)
        @test isnothing(display(t))
        @test isnothing(Trixi.reset_data_structures!(t))
    end
end

@testitem "Unit: TreeMesh" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    @testset "constructors" begin
        mesh = @inferred TreeMesh{1, Trixi.SerialTree{1, Float64}, Float64}(1, 5.0, 2.0,
                                                                            true)
        @test mesh isa TreeMesh

        # Invalid domain length check (TreeMesh expects a hypercube)
        # 2D
        @test_throws ArgumentError TreeMesh((-0.5, 0.0), (1.0, 2.0),
                                            initial_refinement_level = 2,
                                            periodicity = true)
        # 3D
        @test_throws ArgumentError TreeMesh((-0.5, 0.0, -0.2), (1.0, 2.0, 1.5),
                                            initial_refinement_level = 2,
                                            periodicity = true)

        # Keyword-only constructor
        mesh_ref = TreeMesh((-1.0, -1.0), (1.0, 1.0);
                            initial_refinement_level = 2)
        mesh_kw = TreeMesh(; coordinates_min = (-1.0, -1.0),
                           coordinates_max = (1.0, 1.0),
                           refinement_level = 2)
        @test Trixi.ncells(mesh_kw) == Trixi.ncells(mesh_ref)
        @test_throws ArgumentError TreeMesh(; coordinates_min = (-1.0, -1.0),
                                            coordinates_max = (1.0, 1.0, 1.0),
                                            refinement_level = 2)
    end

    @testset "helper functions" begin
        coordinates_min = (-0.5, -0.5, -0.5)
        coordinates_max = (0.5, 0.5, 0.5)

        for ndims in 1:3
            coords_min = coordinates_min[1:ndims]
            coords_max = coordinates_max[1:ndims]
            for ref_level in 0:2
                mesh = TreeMesh(coords_min, coords_max,
                                initial_refinement_level = ref_level,
                                periodicity = true)

                @test @inferred(Trixi.ndims(mesh)) == ndims
                @test @inferred(Trixi.ncells(mesh)) == (2^ndims)^ref_level
            end
        end
    end
end

@testitem "Unit: TreeMeshParallel" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    @testset "partition!" begin
        @testset "mpi_nranks() = 2" begin
            Trixi.mpi_nranks() = 2
            let
                @test @inferred(Trixi.mpi_nranks()) == 2

                mesh = TreeMesh{2, Trixi.ParallelTree{2, Float64}, Float64}(30,
                                                                            (0.0, 0.0),
                                                                            1.0,
                                                                            true)
                # Refine twice
                Trixi.refine!(mesh.tree)
                Trixi.refine!(mesh.tree)

                # allow_coarsening = true
                Trixi.partition!(mesh)
                # Use parent for OffsetArray
                @test parent(mesh.n_cells_by_rank) == [11, 10]
                @test mesh.tree.mpi_ranks[1:21] ==
                      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
                @test parent(mesh.first_cell_by_rank) == [1, 12]

                # allow_coarsening = false
                Trixi.partition!(mesh; allow_coarsening = false)
                @test parent(mesh.n_cells_by_rank) == [11, 10]
                @test mesh.tree.mpi_ranks[1:21] ==
                      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
                @test parent(mesh.first_cell_by_rank) == [1, 12]
            end
            Trixi.mpi_nranks() = Trixi.MPI_SIZE[] # restore the original behavior
        end

        @testset "mpi_nranks() = 3" begin
            Trixi.mpi_nranks() = 3
            let
                @test Trixi.mpi_nranks() == 3

                mesh = TreeMesh{2, Trixi.ParallelTree{2, Float64}, Float64}(100,
                                                                            (0.0, 0.0),
                                                                            1.0,
                                                                            true)
                # Refine twice
                Trixi.refine!(mesh.tree)
                Trixi.refine!(mesh.tree)

                # allow_coarsening = true
                Trixi.partition!(mesh)
                # Use parent for OffsetArray
                @test parent(mesh.n_cells_by_rank) == [11, 5, 5]
                @test mesh.tree.mpi_ranks[1:21] ==
                      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2]
                @test parent(mesh.first_cell_by_rank) == [1, 12, 17]

                # allow_coarsening = false
                Trixi.partition!(mesh; allow_coarsening = false)
                @test parent(mesh.n_cells_by_rank) == [9, 6, 6]
                @test mesh.tree.mpi_ranks[1:21] ==
                      [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
                @test parent(mesh.first_cell_by_rank) == [1, 10, 16]
            end
            Trixi.mpi_nranks() = Trixi.MPI_SIZE[] # restore the original behavior
        end

        @testset "mpi_nranks() = 9" begin
            Trixi.mpi_nranks() = 9
            let
                @test Trixi.mpi_nranks() == 9

                mesh = TreeMesh{2, Trixi.ParallelTree{2, Float64}, Float64}(1000,
                                                                            (0.0, 0.0),
                                                                            1.0,
                                                                            true)
                # Refine twice
                Trixi.refine!(mesh.tree)
                Trixi.refine!(mesh.tree)
                Trixi.refine!(mesh.tree)
                Trixi.refine!(mesh.tree)

                # allow_coarsening = true
                Trixi.partition!(mesh)
                # Use parent for OffsetArray
                @test parent(mesh.n_cells_by_rank) ==
                      [44, 37, 38, 37, 37, 37, 38, 37, 36]
                @test parent(mesh.first_cell_by_rank) ==
                      [1, 45, 82, 120, 157, 194, 231, 269, 306]
            end
            Trixi.mpi_nranks() = Trixi.MPI_SIZE[] # restore the original behavior
        end

        @testset "mpi_nranks() = 3 non-uniform" begin
            Trixi.mpi_nranks() = 3
            let
                @test Trixi.mpi_nranks() == 3

                mesh = TreeMesh{2, Trixi.ParallelTree{2, Float64}, Float64}(100,
                                                                            (0.0, 0.0),
                                                                            1.0,
                                                                            true)
                # Refine whole tree
                Trixi.refine!(mesh.tree)
                # Refine left leaf
                Trixi.refine!(mesh.tree, [2])

                # allow_coarsening = true
                Trixi.partition!(mesh)
                # Use parent for OffsetArray
                @test parent(mesh.n_cells_by_rank) == [6, 1, 2]
                @test mesh.tree.mpi_ranks[1:9] == [0, 0, 0, 0, 0, 0, 1, 2, 2]
                @test parent(mesh.first_cell_by_rank) == [1, 7, 8]

                # allow_coarsening = false
                Trixi.partition!(mesh; allow_coarsening = false)
                @test parent(mesh.n_cells_by_rank) == [5, 2, 2]
                @test mesh.tree.mpi_ranks[1:9] == [0, 0, 0, 0, 0, 1, 1, 2, 2]
                @test parent(mesh.first_cell_by_rank) == [1, 6, 8]
            end
            Trixi.mpi_nranks() = Trixi.MPI_SIZE[] # restore the original behavior
        end

        @testset "not enough ranks" begin
            Trixi.mpi_nranks() = 3
            let
                @test Trixi.mpi_nranks() == 3

                mesh = TreeMesh{2, Trixi.ParallelTree{2, Float64}, Float64}(100,
                                                                            (0.0, 0.0),
                                                                            1.0,
                                                                            true)

                # Only one leaf
                @test_throws AssertionError("Too many ranks to properly partition the mesh!") Trixi.partition!(mesh)

                # Refine to 4 leaves
                Trixi.refine!(mesh.tree)

                # All four leaves will need to be on one rank to allow coarsening
                @test_throws AssertionError("Too many ranks to properly partition the mesh!") Trixi.partition!(mesh)
                @test_nowarn Trixi.partition!(mesh; allow_coarsening = false)
            end
            Trixi.mpi_nranks() = Trixi.MPI_SIZE[] # restore the original behavior
        end
    end
end

@testitem "Unit: curved mesh" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    @testset "calc_jacobian_matrix" begin
        @testset "identity map" begin
            basis = LobattoLegendreBasis(5)
            nodes = Trixi.get_nodes(basis)
            jacobian_matrix = Array{Float64, 5}(undef, 2, 2, 6, 6, 1)

            node_coordinates = Array{Float64, 4}(undef, 2, 6, 6, 1)
            node_coordinates[1, :, :, 1] .= [nodes[i] for i in 1:6, j in 1:6]
            node_coordinates[2, :, :, 1] .= [nodes[j] for i in 1:6, j in 1:6]
            expected = zeros(2, 2, 6, 6, 1)
            expected[1, 1, :, :, 1] .= 1
            expected[2, 2, :, :, 1] .= 1
            @test Trixi.calc_jacobian_matrix!(jacobian_matrix, 1, node_coordinates,
                                              basis) ≈ expected
        end

        @testset "maximum exact polydeg" begin
            basis = LobattoLegendreBasis(3)
            nodes = Trixi.get_nodes(basis)
            jacobian_matrix = Array{Float64, 5}(undef, 2, 2, 4, 4, 1)

            # f(x, y) = [x^3, xy^2]
            node_coordinates = Array{Float64, 4}(undef, 2, 4, 4, 1)
            node_coordinates[1, :, :, 1] .= [nodes[i]^3 for i in 1:4, j in 1:4]
            node_coordinates[2, :, :, 1] .= [nodes[i] * nodes[j]^2
                                             for i in 1:4, j in 1:4]

            # Df(x, y) = [3x^2 0;
            #              y^2 2xy]
            expected = zeros(2, 2, 4, 4, 1)
            expected[1, 1, :, :, 1] .= [3 * nodes[i]^2 for i in 1:4, j in 1:4]
            expected[2, 1, :, :, 1] .= [nodes[j]^2 for i in 1:4, j in 1:4]
            expected[2, 2, :, :, 1] .= [2 * nodes[i] * nodes[j] for i in 1:4, j in 1:4]
            @test Trixi.calc_jacobian_matrix!(jacobian_matrix, 1, node_coordinates,
                                              basis) ≈ expected
        end
    end
end

@testitem "Unit: interpolation" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    @testset "nodes and weights" begin
        @test Trixi.gauss_nodes_weights(1) == ([0.0], [2.0])

        @test Trixi.gauss_nodes_weights(2)[1] ≈ [-1 / sqrt(3), 1 / sqrt(3)]
        @test Trixi.gauss_nodes_weights(2)[2] == [1.0, 1.0]

        @test Trixi.gauss_nodes_weights(3)[1] ≈ [-sqrt(3 / 5), 0.0, sqrt(3 / 5)]
        @test Trixi.gauss_nodes_weights(3)[2] ≈ [5 / 9, 8 / 9, 5 / 9]
    end

    @testset "boundary interpolation" begin
        for p in 1:7
            basis = LobattoLegendreBasis(p)
            nodes = basis.nodes
            weights = basis.weights

            L_minus1 = Trixi.calc_L(-1.0, nodes, weights)
            Lhat_minus1 = Trixi.calc_Lhat(L_minus1, weights)
            @test basis.inverse_weights[1] == Lhat_minus1[1]

            L_plus1 = Trixi.calc_L(1.0, nodes, weights)
            Lhat_plus1 = Trixi.calc_Lhat(L_plus1, weights)
            @test basis.inverse_weights[p + 1] == Lhat_plus1[p + 1]
        end
    end

    @testset "multiply_dimensionwise" begin
        nodes_in = [0.0, 0.5, 1.0]
        nodes_out = [0.0, 1 / 3, 2 / 3, 1.0]
        matrix = Trixi.polynomial_interpolation_matrix(nodes_in, nodes_out)
        data_in = [3.0 4.5 6.0]
        @test isapprox(Trixi.multiply_dimensionwise(matrix, data_in), [3.0 4.0 5.0 6.0])

        n_vars = 3
        size_in = 2
        size_out = 3
        matrix = randn(size_out, size_in)
        # 1D
        data_in = randn(n_vars, size_in)
        data_out = Trixi.multiply_dimensionwise_naive(matrix, data_in)
        @test isapprox(data_out, Trixi.multiply_dimensionwise(matrix, data_in))
        # 2D
        data_in = randn(n_vars, size_in, size_in)
        data_out = Trixi.multiply_dimensionwise_naive(matrix, data_in)
        @test isapprox(data_out, Trixi.multiply_dimensionwise(matrix, data_in))
        # 3D
        data_in = randn(n_vars, size_in, size_in, size_in)
        data_out = Trixi.multiply_dimensionwise_naive(matrix, data_in)
        @test isapprox(data_out, Trixi.multiply_dimensionwise(matrix, data_in))
    end
end

@testitem "Unit: L2 projection" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    @testset "calc_reverse_upper for LGL" begin
        @test isapprox(Trixi.calc_reverse_upper(2, Val(:gauss_lobatto)),
                       [[0.25, 0.25] [0.0, 0.5]])
    end
    @testset "calc_reverse_lower for LGL" begin
        @test isapprox(Trixi.calc_reverse_lower(2, Val(:gauss_lobatto)),
                       [[0.5, 0.0] [0.25, 0.25]])
    end
end

@testitem "Unit: GaussLegendreBasis" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    basis = GaussLegendreBasis(3)
    @test nnodes(basis) == 4
    @test_nowarn show(stdout, "text/plain", basis)

    solution_analyzer = Trixi.SolutionAnalyzer(basis)
    @test nnodes(solution_analyzer) == 7
    @test_nowarn show(stdout, "text/plain", solution_analyzer)
end

@testitem "Unit: Positivity limiter for AMRCallback" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Initial condition with simple discontinuity
    @inline function initial_condition_discontinuity(x, t, equations)
        if x[1] < -0.5
            scalar = 1.0e-4
        else
            scalar = 1.0
        end

        return SVector(scalar)
    end

    # Set up variables used for 1D, 2D and 3D tests
    solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs)

    adaptor = Trixi.AdaptorL2(solver.basis)
    limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6,),
                                                   variables = (first,))

    @testset "1D" begin
        equations = LinearScalarAdvectionEquation1D(1.0)

        coordinates_min = (-1.0,)
        coordinates_max = (1.0,)
        mesh = TreeMesh(coordinates_min, coordinates_max,
                        initial_refinement_level = 0,
                        periodicity = true)

        semi = SemidiscretizationHyperbolic(mesh, equations,
                                            initial_condition_discontinuity, solver,
                                            boundary_conditions = boundary_condition_periodic)

        u_ode = compute_coefficients(initial_condition_discontinuity, 0.0, semi)

        # Refinement
        elements_to_refine = [1]
        Trixi.refine!(mesh.tree, elements_to_refine)

        Trixi.refine!(u_ode, adaptor, mesh, equations, solver, semi.cache,
                      elements_to_refine, limiter!)

        @test all(u_ode .>= 0.0)

        # Coarsening
        u_ode = compute_coefficients(initial_condition_discontinuity, 0.0, semi)

        parents_to_coarsen = [1]
        Trixi.coarsen!(mesh.tree, parents_to_coarsen)

        elements_to_coarsen = collect(1:2)
        Trixi.coarsen!(u_ode, adaptor, mesh, equations, solver, semi.cache,
                       elements_to_coarsen, limiter!)

        @test all(u_ode .>= 0.0)
    end

    @testset "2D" begin
        equations = LinearScalarAdvectionEquation2D((0.2, -0.7))

        coordinates_min = (-1.0, -1.0)
        coordinates_max = (1.0, 1.0)
        mesh = TreeMesh(coordinates_min, coordinates_max,
                        initial_refinement_level = 0,
                        periodicity = true)

        semi = SemidiscretizationHyperbolic(mesh, equations,
                                            initial_condition_discontinuity, solver,
                                            boundary_conditions = boundary_condition_periodic)

        u_ode = compute_coefficients(initial_condition_discontinuity, 0.0, semi)

        # Refinement
        elements_to_refine = [1]
        Trixi.refine!(mesh.tree, elements_to_refine)

        Trixi.refine!(u_ode, adaptor, mesh, equations, solver, semi.cache,
                      elements_to_refine, limiter!)

        @test all(u_ode .>= 0.0)

        # Coarsening
        u_ode = compute_coefficients(initial_condition_discontinuity, 0.0, semi)

        parents_to_coarsen = [1]
        Trixi.coarsen!(mesh.tree, parents_to_coarsen)

        elements_to_coarsen = collect(1:4)
        Trixi.coarsen!(u_ode, adaptor, mesh, equations, solver, semi.cache,
                       elements_to_coarsen, limiter!)

        @test all(u_ode .>= 0.0)
    end

    @testset "3D" begin
        equations = LinearScalarAdvectionEquation3D((0.2, -0.7, 0.5))

        coordinates_min = (-1.0, -1.0, -1.0)
        coordinates_max = (1.0, 1.0, 1.0)
        mesh = TreeMesh(coordinates_min, coordinates_max,
                        initial_refinement_level = 0,
                        periodicity = true)

        semi = SemidiscretizationHyperbolic(mesh, equations,
                                            initial_condition_discontinuity, solver,
                                            boundary_conditions = boundary_condition_periodic)

        u_ode = compute_coefficients(initial_condition_discontinuity, 0.0, semi)

        # Refinement
        elements_to_refine = [1]
        Trixi.refine!(mesh.tree, elements_to_refine)

        Trixi.refine!(u_ode, adaptor, mesh, equations, solver, semi.cache,
                      elements_to_refine, limiter!)

        @test all(u_ode .>= 0.0)

        # Coarsening
        u_ode = compute_coefficients(initial_condition_discontinuity, 0.0, semi)

        parents_to_coarsen = [1]
        Trixi.coarsen!(mesh.tree, parents_to_coarsen)

        elements_to_coarsen = collect(1:8)
        Trixi.coarsen!(u_ode, adaptor, mesh, equations, solver, semi.cache,
                       elements_to_coarsen, limiter!)

        @test all(u_ode .>= 0.0)
    end
end

@testitem "Unit: containers" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Set up mock container
    mutable struct MyContainer <: Trixi.AbstractContainer
        data::Vector{Int}
        capacity::Int
        length::Int
        dummy::Int
    end
    function MyContainer(data, capacity)
        c = MyContainer(Vector{Int}(undef, capacity + 1), capacity, length(data),
                        capacity + 1)
        c.data[eachindex(data)] .= data
        return c
    end
    MyContainer(data::AbstractArray) = MyContainer(data, length(data))
    Trixi.invalidate!(c::MyContainer, first, last) = (c.data[first:last] .= 0; c)
    function Trixi.raw_copy!(target::MyContainer, source::MyContainer, first, last,
                             destination)
        Trixi.copy_data!(target.data, source.data, first, last, destination)
        return target
    end
    Trixi.move_connectivity!(c::MyContainer, first, last, destination) = c
    Trixi.delete_connectivity!(c::MyContainer, first, last) = c
    function Trixi.reset_data_structures!(c::MyContainer)
        c.data = Vector{Int}(undef, c.capacity + 1)
        return c
    end
    function Base.:(==)(c1::MyContainer, c2::MyContainer)
        return (c1.capacity == c2.capacity &&
                c1.length == c2.length &&
                c1.dummy == c2.dummy &&
                c1.data[1:(c1.length)] == c2.data[1:(c2.length)])
    end

    @testset "size" begin
        c = MyContainer([1, 2, 3])
        @test size(c) == (3,)
    end

    @testset "resize!" begin
        c = MyContainer([1, 2, 3])
        @test length(resize!(c, 2)) == 2
    end

    @testset "copy!" begin
        c1 = MyContainer([1, 2, 3])
        c2 = MyContainer([4, 5])
        @test Trixi.copy!(c1, c2, 2, 1, 2) == MyContainer([1, 2, 3]) # no-op

        c1 = MyContainer([1, 2, 3])
        c2 = MyContainer([4, 5])
        @test Trixi.copy!(c1, c2, 1, 2, 2) == MyContainer([1, 4, 5])

        c1 = MyContainer([1, 2, 3])
        @test Trixi.copy!(c1, c2, 1, 2) == MyContainer([1, 4, 3])

        c1 = MyContainer([1, 2, 3])
        @test Trixi.copy!(c1, 2, 3, 1) == MyContainer([2, 3, 3])

        c1 = MyContainer([1, 2, 3])
        @test Trixi.copy!(c1, 1, 3) == MyContainer([1, 2, 1])
    end

    @testset "move!" begin
        c = MyContainer([1, 2, 3])
        @test Trixi.move!(c, 1, 1) == MyContainer([1, 2, 3]) # no-op

        c = MyContainer([1, 2, 3])
        @test Trixi.move!(c, 1, 2) == MyContainer([0, 1, 3])
    end

    @testset "remove_shift!" begin
        c = MyContainer([1, 2, 3, 4])
        @test Trixi.remove_shift!(c, 2, 1) == MyContainer([1, 2, 3, 4]) # no-op

        c = MyContainer([1, 2, 3, 4])
        @test Trixi.remove_shift!(c, 2, 2) == MyContainer([1, 3, 4], 4)

        c = MyContainer([1, 2, 3, 4])
        @test Trixi.remove_shift!(c, 2) == MyContainer([1, 3, 4], 4)
    end
end

@testitem "Unit: example elixirs" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    @test basename(examples_dir()) == "examples"
    @test !isempty(get_examples())
    @test endswith(default_example(), "elixir_advection_basic.jl")
end

@testitem "Unit: HLL flux with vanishing wave speed estimates (#502)" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    equations = CompressibleEulerEquations1D(1.4)
    u = SVector(1.0, 0.0, 0.0)
    @test !any(isnan, flux_hll(u, u, 1, equations))
end

@testitem "Unit: DG L2 mortar container debug output" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    c2d = Trixi.TreeL2MortarContainer2D{Float64}(1, 1, 1)
    @test isnothing(display(c2d))
    c3d = Trixi.TreeL2MortarContainer3D{Float64}(1, 1, 1)
    @test isnothing(display(c3d))
end

@testitem "Unit: TreeContainer1D nnodes(container)" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    capacity = 42
    n_variables = 9

    interface_container = Trixi.TreeInterfaceContainer1D{Float64}(capacity, n_variables)
    @test nnodes(interface_container) == 1

    boundary_container = Trixi.TreeBoundaryContainer1D{Float64, Float64}(capacity,
                                                                         n_variables)
    @test nnodes(boundary_container) == 1
end

@testitem "Unit: Printing indicators/controllers" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Note: Constructing indicators/controllers using the parameters below doesn't make sense. It's
    # just useful to run basic tests of `show` methods.

    c = ControllerThreeLevelCombined(1, 2, 3, 10.0, 11.0, 12.0, "primary", "secondary",
                                     "cache")
    @test_nowarn show(stdout, c)

    indicator_hg = IndicatorHennemannGassner(1.0, 0.0, true, "variable", "cache")
    @test_nowarn show(stdout, indicator_hg)

    limiter_idp = SubcellLimiterIDP(true, [1], true, [1], ["variable"], 0.1,
                                    true, [(entropy_guermond_etal, min)], "cache",
                                    1, (1.0, 1.0), 1.0)
    @test_nowarn show(stdout, limiter_idp)

    indicator_loehner = IndicatorLöhner(1.0, "variable", (; cache = nothing))
    @test_nowarn show(stdout, indicator_loehner)

    indicator_max = IndicatorMax("variable", (; cache = nothing))
    @test_nowarn show(stdout, indicator_max)

    indicator_ec = IndicatorEntropyCorrection(CompressibleEulerEquations1D(1.4),
                                              LobattoLegendreBasis(3))
    @test_nowarn show(stdout, indicator_ec)

    # test Base.show for PositivityPreservingLimiterLiuZhang
    equations = LinearScalarAdvectionEquation1D(1.0)
    solver = DGSEM(polydeg = 3)
    mesh = TreeMesh(-1.0, 1.0, initial_refinement_level = 1, periodicity = true)
    semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_constant,
                                        solver;
                                        boundary_conditions = boundary_condition_periodic)
    local_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (1e-3,),
                                                         variables = (first,))
    global_limiter! = PositivityPreservingLimiterLiuZhang(local_limiter!,
                                                          semi;
                                                          record_davis_yin_iterations = true)
    @test_nowarn show(stdout, global_limiter!)
    @test_nowarn show(stdout, "text/plain", global_limiter!)
    @test_nowarn show(IOContext(IOBuffer(), :compact => true), MIME"text/plain"(),
                      global_limiter!)
    @test_nowarn show(IOContext(IOBuffer(), :compact => false), MIME"text/plain"(),
                      global_limiter!)
end

@testitem "Unit: LBM 2D constructor" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Neither Mach number nor velocity set
    @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = nothing, Re = 1000)
    # Both Mach number and velocity set
    @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000,
                                                            u0 = 1.0)
    # Neither Reynolds number nor viscosity set
    @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = nothing)
    # Both Reynolds number and viscosity set
    @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000,
                                                            nu = 1.0)

    # No non-dimensional values set
    @test LatticeBoltzmannEquations2D(Ma = nothing, Re = nothing, u0 = 1.0,
                                      nu = 1.0) isa
          LatticeBoltzmannEquations2D
end

@testitem "Unit: LBM 3D constructor" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Neither Mach number nor velocity set
    @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = nothing, Re = 1000)
    # Both Mach number and velocity set
    @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000,
                                                            u0 = 1.0)
    # Neither Reynolds number nor viscosity set
    @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = nothing)
    # Both Reynolds number and viscosity set
    @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000,
                                                            nu = 1.0)

    # No non-dimensional values set
    @test LatticeBoltzmannEquations3D(Ma = nothing, Re = nothing, u0 = 1.0,
                                      nu = 1.0) isa
          LatticeBoltzmannEquations3D
end

@testitem "Unit: LBM 2D functions" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Set up LBM struct and dummy distribution
    equation = LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000)
    u = Trixi.equilibrium_distribution(1, 2, 3, equation)

    # Component-wise velocity
    @test isapprox(Trixi.velocity(u, 1, equation), 2)
    @test isapprox(Trixi.velocity(u, 2, equation), 3)
end

@testitem "Unit: LBM 3D functions" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Set up LBM struct and dummy distribution
    equation = LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000)
    u = Trixi.equilibrium_distribution(1, 2, 3, 4, equation)

    # Component-wise velocity
    @test isapprox(velocity(u, 1, equation), 2)
    @test isapprox(velocity(u, 2, equation), 3)
    @test isapprox(velocity(u, 3, equation), 4)
end

@testitem "Unit: LBMCollisionCallback" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Printing of LBM collision callback
    callback = LBMCollisionCallback()
    @test_nowarn show(stdout, callback)
    println()
    @test_nowarn show(stdout, "text/plain", callback)
    println()
end

@testitem "Unit: Acoustic perturbation 2D varnames" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    v_mean_global = (0.0, 0.0)
    c_mean_global = 1.0
    rho_mean_global = 1.0
    equations = AcousticPerturbationEquations2D(v_mean_global, c_mean_global,
                                                rho_mean_global)

    @test Trixi.varnames(cons2state, equations) ==
          ("v1_prime", "v2_prime", "p_prime_scaled")
    @test Trixi.varnames(cons2mean, equations) ==
          ("v1_mean", "v2_mean", "c_mean", "rho_mean")
end

@testitem "Unit: Euler conversion between conservative/entropy variables" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    rho, v1, v2, v3, p = 1.0, 0.1, 0.2, 0.3, 2.0

    let equations = CompressibleEulerEquations1D(1.4)
        cons_vars = prim2cons(SVector(rho, v1, p), equations)
        entropy_vars = cons2entropy(cons_vars, equations)
        @test cons_vars ≈ entropy2cons(entropy_vars, equations)

        # test tuple args
        cons_vars = prim2cons((rho, v1, p), equations)
        entropy_vars = cons2entropy(cons_vars, equations)
        @test cons_vars ≈ entropy2cons(entropy_vars, equations)
    end

    # Test PassiveTracerEquations
    let flow_equations = CompressibleEulerEquations1D(1.4)
        equations = PassiveTracerEquations(flow_equations, n_tracers = 2)
        xi1, xi2 = 0.4, 0.5
        cons_ref = SVector(rho, rho * v1, p / 0.4 + 0.5 * (rho * v1 * v1), rho * xi1,
                           rho * xi2)
        cons_test = prim2cons(SVector(rho, v1, p, xi1, xi2), equations)
        @test cons_test ≈ cons_ref
        prim_test = cons2prim(cons_test, equations)
        @test prim_test ≈ SVector(rho, v1, p, xi1, xi2)
        flow_entropy = cons2entropy(cons_ref, flow_equations)

        entropy_ref = SVector(flow_entropy[1] - (xi1^2 + xi2^2),
                              (flow_entropy[i] for i in 2:nvariables(flow_equations))...,
                              2 * xi1, 2 * xi2)
        entropy_test = cons2entropy(cons_test, equations)
        @test entropy_test ≈ entropy_ref

        # Also test density, pressure, density_pressure and entropy here because there is currently
        # no specific space for testing them (e.g., in the other equations)
        @test density(cons_test, equations) ≈ rho
        @test pressure(cons_test, equations) ≈ p
        @test density_pressure(cons_test, equations) ≈ rho * p
        @test entropy(cons_test, equations) ≈
              entropy(cons_ref, flow_equations) + rho * (xi1^2 + xi2^2)

        tracers_ = Trixi.tracers(cons_test, equations)
        @test tracers_ ≈ SVector(xi1, xi2)
        rho_tracers_ = Trixi.rho_tracers(cons_test, equations)
        @test rho_tracers_ ≈ SVector(rho * xi1, rho * xi2)
    end

    let equations = CompressibleEulerEquations2D(1.4)
        cons_vars = prim2cons(SVector(rho, v1, v2, p), equations)
        entropy_vars = cons2entropy(cons_vars, equations)
        @test cons_vars ≈ entropy2cons(entropy_vars, equations)

        # test tuple args
        cons_vars = prim2cons((rho, v1, v2, p), equations)
        entropy_vars = cons2entropy(cons_vars, equations)
        @test cons_vars ≈ entropy2cons(entropy_vars, equations)
    end

    let equations = CompressibleEulerEquations3D(1.4)
        cons_vars = prim2cons(SVector(rho, v1, v2, v3, p), equations)
        entropy_vars = cons2entropy(cons_vars, equations)
        @test cons_vars ≈ entropy2cons(entropy_vars, equations)

        # test tuple args
        cons_vars = prim2cons((rho, v1, v2, v3, p), equations)
        entropy_vars = cons2entropy(cons_vars, equations)
        @test cons_vars ≈ entropy2cons(entropy_vars, equations)
    end
end

@testitem "Unit: Navier-Stokes conversion between conservative/primitive variables" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    rho, v1, v2, v3, p = 2.0, 0.1, 0.2, 0.3, 4.0
    mu, Prandtl = 0.01, 0.72

    let equations_hyperbolic = CompressibleEulerEquations1D(1.4)
        equations = CompressibleNavierStokesDiffusion1D(equations_hyperbolic;
                                                        mu, Prandtl)
        prim_vars = SVector(rho, v1, p)
        cons_vars = prim2cons(prim_vars, equations)
        @test prim_vars ≈ cons2prim(cons_vars, equations)
        @test cons_vars ≈ prim2cons(cons2prim(cons_vars, equations), equations)
    end

    let equations_hyperbolic = CompressibleEulerEquations2D(1.4)
        equations = CompressibleNavierStokesDiffusion2D(equations_hyperbolic;
                                                        mu, Prandtl)
        prim_vars = SVector(rho, v1, v2, p)
        cons_vars = prim2cons(prim_vars, equations)
        @test prim_vars ≈ cons2prim(cons_vars, equations)
        @test cons_vars ≈ prim2cons(cons2prim(cons_vars, equations), equations)
    end

    let equations_hyperbolic = CompressibleEulerEquations3D(1.4)
        equations = CompressibleNavierStokesDiffusion3D(equations_hyperbolic;
                                                        mu, Prandtl)
        prim_vars = SVector(rho, v1, v2, v3, p)
        cons_vars = prim2cons(prim_vars, equations)
        @test prim_vars ≈ cons2prim(cons_vars, equations)
        @test cons_vars ≈ prim2cons(cons2prim(cons_vars, equations), equations)
    end
end

@testitem "Unit: LaplaceDiffusionEntropyVariables apply_jacobian_entropy2cons" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    rho, v1, v2, v3, p = 1.0, 0.1, -0.2, 0.3, 2.0

    for (equations_parabolic, prim, dw) in ((LaplaceDiffusionEntropyVariables1D(0.01,
                                                                                CompressibleEulerEquations1D(1.4)),
                                             SVector(rho, v1, p),
                                             SVector(1.0, 0.1, 2.0)),
                                            (LaplaceDiffusionEntropyVariables2D(0.01,
                                                                                CompressibleEulerEquations2D(1.4)),
                                             SVector(rho, v1, v2, p),
                                             SVector(1.0, 0.1, -0.2, 2.0)),
                                            (LaplaceDiffusionEntropyVariables3D(0.01,
                                                                                CompressibleEulerEquations3D(1.4)),
                                             SVector(rho, v1, v2, v3, p),
                                             SVector(1.0, 0.1, -0.2, 0.3, 2.0)))
        equations = equations_parabolic.equations_hyperbolic
        w = cons2entropy(prim2cons(prim, equations), equations)
        jvp_specialized = Trixi.apply_jacobian_entropy2cons(dw, w, equations_parabolic)
        jvp_ad = invoke(Trixi.apply_jacobian_entropy2cons,
                        Tuple{typeof(dw), typeof(w), Trixi.AbstractEquations},
                        dw, w, equations)
        @test jvp_specialized ≈ jvp_ad
    end
end

# It is for many equations possible to compute ρ ⋅ p more efficiently
# than computing the pressure (and density if needed) separately and then multiplying.
# This is due to the computation of the kinetic energy term, which usually involves
# dividing the squared momenta by the density, an operation that can be avoided
# when computing the product ρ ⋅ p directly.
@testitem "Unit: Test density_pressure" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    let equations = CompressibleEulerEquations1D(5 / 3)
        u = initial_condition_density_wave(SVector(1.0), 3.0, equations)
        rho = density(u, equations)
        p = pressure(u, equations)
        rho_p = density_pressure(u, equations)
        @test rho * p ≈ rho_p
    end

    let equations = CompressibleEulerEquations2D(2.0)
        u = initial_condition_eoc_test_coupled_euler_gravity(SVector(0.666, 0.25), 0.1,
                                                             equations)
        rho = density(u, equations)
        p = pressure(u, equations)
        rho_p = density_pressure(u, equations)
        @test rho * p ≈ rho_p
    end

    let equations = CompressibleEulerEquations3D(1.4)
        u = initial_condition_convergence_test(SVector(0.5, 0.1, -0.2), 1.5, equations)
        rho = density(u, equations)
        p = pressure(u, equations)
        rho_p = density_pressure(u, equations)
        @test rho * p ≈ rho_p
    end

    let equations = CompressibleEulerEquationsQuasi1D(1.4)
        u = initial_condition_convergence_test(SVector(2.0), 5.0, equations)
        rho = density(u, equations)
        p = pressure(u, equations)
        rho_p = density_pressure(u, equations)
        @test rho * p ≈ rho_p
    end

    let equations = IdealGlmMhdEquations1D(5 / 3)
        u = initial_condition_convergence_test(SVector(-1.0), 7.0, equations)
        rho = density(u, equations)
        p = pressure(u, equations)
        rho_p = density_pressure(u, equations)
        @test rho * p ≈ rho_p
    end

    let equations = IdealGlmMhdEquations2D(5 / 3)
        u = initial_condition_convergence_test(SVector(-1.0, 0.5), 0.1, equations)
        rho = density(u, equations)
        p = pressure(u, equations)
        rho_p = density_pressure(u, equations)
        @test rho * p ≈ rho_p
    end

    let equations = IdealGlmMhdEquations3D(5 / 3)
        u = initial_condition_convergence_test(SVector(-1.0, 0.5, 0.2), 0.8, equations)
        rho = density(u, equations)
        p = pressure(u, equations)
        rho_p = density_pressure(u, equations)
        @test rho * p ≈ rho_p
    end

    let equations = CompressibleEulerMulticomponentEquations1D(gammas = (1.4, 1.4),
                                                               gas_constants = (0.4,
                                                                                0.4))
        u = initial_condition_convergence_test(SVector(1.0), 42.0, equations)
        rho = density(u, equations)
        p = pressure(u, equations)
        rho_p = density_pressure(u, equations)
        @test rho * p ≈ rho_p
    end

    let equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648),
                                                               gas_constants = (0.287,
                                                                                1.578))
        u = initial_condition_convergence_test(SVector(1.0, 0.1), 0.42, equations)
        rho = density(u, equations)
        p = pressure(u, equations)
        rho_p = density_pressure(u, equations)
        @test rho * p ≈ rho_p
    end

    let equations = IdealGlmMhdMulticomponentEquations1D(gammas = (2.0, 2.0, 2.0),
                                                         gas_constants = (2.0, 2.0,
                                                                          2.0))
        u = initial_condition_weak_blast_wave(SVector(0.5, 0.1), 0.0, equations)
        rho = density(u, equations)
        p = pressure(u, equations)
        rho_p = density_pressure(u, equations)
        @test rho * p ≈ rho_p
    end

    let equations = IdealGlmMhdMulticomponentEquations2D(gammas = (5 / 3, 5 / 3, 5 / 3),
                                                         gas_constants = (2.08, 2.08,
                                                                          2.08))
        u = initial_condition_convergence_test(SVector(-0.5, 0.1), 0.666, equations)
        rho = density(u, equations)
        p = pressure(u, equations)
        rho_p = density_pressure(u, equations)
        @test rho * p ≈ rho_p
    end
end

@testitem "Unit: hardened boundary_condition_slip_wall" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    let equations = CompressibleEulerEquations1D(1.4)
        # this state results in base < 0 in the implementation of
        # boundary_condition_slip_wall.
        u_inner = prim2cons(SVector(1.4, -6, 1), equations)

        x, t = 0, 0
        orientation = 1
        direction = 2 # even direction: v_normal stays <= 0
        surface_flux_function = flux_lax_friedrichs
        flux = boundary_condition_slip_wall(u_inner, orientation, direction, x, t,
                                            surface_flux_function, equations)

        # boundary_condition_slip_wall should return exactly zero
        @test flux == zero(flux)
    end

    let equations = CompressibleEulerEquations2D(1.4)
        # this state results in base < 0 in the implementation of
        # boundary_condition_slip_wall.
        u_inner = prim2cons(SVector(1.4, -6, 0, 1), equations)
        normal_direction = SVector(1.0, 0.0)

        x, t = 0, 0
        surface_flux_function = flux_lax_friedrichs
        flux = boundary_condition_slip_wall(u_inner, normal_direction, x, t,
                                            surface_flux_function, equations)

        # boundary_condition_slip_wall should return exactly zero
        @test flux == zero(flux)
    end

    let equations = CompressibleEulerEquations3D(1.4)
        # this state results in base < 0 in the implementation of
        # boundary_condition_slip_wall.
        u_inner = prim2cons(SVector(1.4, -6, 0, 0, 1), equations)
        normal_direction = SVector(1.0, 0.0, 0.0)

        x, t = 0, 0
        surface_flux_function = flux_lax_friedrichs
        flux = boundary_condition_slip_wall(u_inner, normal_direction, x, t,
                                            surface_flux_function, equations)

        # boundary_condition_slip_wall should return exactly zero
        @test flux == zero(flux)
    end
end

@testitem "Unit: Helmholtz ideal gas equation of state (AD vs analytical)" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    # Closed forms for ideal gas from Klein et al., Appendix E (E.1), in (V, T) variables
    function ideal_gas_analytical_helmholtz(V, T, eos::HelmholtzIdealGas)
        alpha = inv(eos.gamma - 1)
        p = eos.R * T / V
        s = eos.R * (1 + alpha + log((T^alpha) * V))
        e = eos.R * T * alpha
        return (; p, s, e)
    end

    rho = 1.225
    T = 300.15
    V = inv(rho)
    eos = HelmholtzIdealGas()
    ref = ideal_gas_analytical_helmholtz(V, T, eos)
    @test isapprox(pressure(V, T, eos), ref.p)
    @test isapprox(Trixi.entropy_specific(V, T, eos), ref.s)
    @test isapprox(energy_internal_specific(V, T, eos), ref.e)

    ig = IdealGas(1.4, eos.R)
    @test Trixi.speed_of_sound(V, T, eos) ≈ Trixi.speed_of_sound(V, T, ig)
    c_direct = Trixi.speed_of_sound(V, T, eos)
    c_fallback = invoke(Trixi.speed_of_sound,
                        Tuple{typeof(V), typeof(T), Trixi.AbstractHelmholtzEOS},
                        V, T, eos)
    @test c_direct ≈ c_fallback
    @test temperature(V, ref.e, eos) ≈ T
    e_h = energy_internal_specific(V, T, eos)
    p_h = pressure(V, T, eos)
    s_h = Trixi.entropy_specific(V, T, eos)
    @test Trixi.gibbs_free_energy(V, T, eos) ≈ e_h + p_h * V - T * s_h
end

@testitem "Unit: boundary_condition_do_nothing" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    rho, v1, v2, p = 1.0, 0.1, 0.2, 0.3, 2.0

    let equations = CompressibleEulerEquations2D(1.4)
        u = prim2cons(SVector(rho, v1, v2, p), equations)
        x = SVector(1.0, 2.0)
        t = 0.5
        surface_flux = flux_lax_friedrichs

        outward_direction = SVector(0.2, -0.3)
        @test flux(u, outward_direction, equations) ≈
              boundary_condition_do_nothing(u, outward_direction, x, t, surface_flux,
                                            equations)

        orientation = 2
        direction = 4
        @test flux(u, orientation, equations) ≈
              boundary_condition_do_nothing(u, orientation, direction, x, t,
                                            surface_flux, equations)
    end
end

@testitem "Unit: boundary_condition_do_nothing_non_conservative" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    rho, v1, v2, v3, p, B1, B2, B3, psi = 1.0, 0.1, 0.2, 0.3, 1.0, 0.0,
                                          40.0 / sqrt(4.0 * pi), 0.0, 0.0

    let equations = IdealGlmMhdEquations2D(1.4, initial_c_h = 1.0)
        u = prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
        x = SVector(1.0, 2.0)
        t = 0.5
        surface_fluxes = (flux_lax_friedrichs, flux_nonconservative_powell)

        outward_direction = SVector(0.2, 0.3)

        @test all(isapprox(x, y)
                  for (x, y) in zip(ntuple(i -> surface_fluxes[i](u, u,
                                                                  outward_direction,
                                                                  equations), 2),
                                    boundary_condition_do_nothing(u, outward_direction,
                                                                  x, t, surface_fluxes,
                                                                  equations)))

        orientation = 2
        direction = 4

        @test all(isapprox(x, y)
                  for (x, y) in zip(ntuple(i -> surface_fluxes[i](u, u, orientation,
                                                                  equations), 2),
                                    boundary_condition_do_nothing(u, orientation,
                                                                  direction, x, t,
                                                                  surface_fluxes,
                                                                  equations)))
    end
end

@testitem "Unit: Reproducing ideal gas with ThermallyPerfectGas9PolyFit" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    R_specific = 287.0509010514002 # [J/(kg*K)]
    p_ref = 100000.0 # [Pa]
    T_ref = 298.15 # [K]

    gamma_target = 1.4
    cp_over_R = gamma_target / (gamma_target - 1)

    temp_bounds = SVector(eps(Float64), typemax(Float64))  # single wide interval, avoid eps/typemax edge cases
    a = Trixi.SMatrix{9, 1}([0.0, 0.0, cp_over_R, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    eos = ThermallyPerfectGas9PolyFit(R_specific = R_specific,
                                      temperature_bounds = temp_bounds,
                                      coefficients = a,
                                      p_ref = p_ref,
                                      T_ref = T_ref)

    rho = 1.255 # [kg/m^3]
    V = 1 / rho

    # 1. cp, cv, gamma match the ideal-gas properties
    cp = Trixi.heat_capacity_constant_pressure(T_ref, eos)
    cv = Trixi.heat_capacity_constant_volume(V, T_ref, eos)
    @test cp ≈ gamma_target / (gamma_target - 1) * R_specific
    @test cv ≈ 1 / (gamma_target - 1) * R_specific
    @test cp / cv ≈ gamma_target

    # 2. pressure matches ideal gas law
    @test pressure(V, T_ref, eos) ≈ rho * R_specific * T_ref

    # 3. internal energy matches u = cv * T (calorically perfect gas, up to a reference offset)
    T_test = 400.0
    e_internal1 = Trixi.energy_internal_specific(V, T_ref, eos)
    e_internal2 = Trixi.energy_internal_specific(V, T_test, eos)
    @test (e_internal2 - e_internal1) ≈ cv * (T_test - T_ref)

    # 4. entropy difference matches ideal-gas relation
    #    Δs = cp*ln(T2/T1) - R_specific*ln(p2/p1)
    p1 = pressure(V, T_ref, eos)
    p2 = pressure(V, T_test, eos)  # same V, different T -> different p
    s1 = Trixi.entropy_specific(V, T_ref, eos)
    s2 = Trixi.entropy_specific(V, T_test, eos)
    @test (s2 - s1) ≈ cp * log(T_test / T_ref) - R_specific * log(p2 / p1)

    # 5. check speed of sound
    @test Trixi.speed_of_sound(V, T_ref, eos) ≈
          sqrt(gamma_target * pressure(V, T_ref, eos) * V)
end

@testitem "Unit: Test consistency (fluxes, entropy/cons2entropy) for NonIdealCompressibleEulerEquations1D" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    eos = VanDerWaals(; a = 10, b = 0.01, R = 287, gamma = 1.4)
    equations = NonIdealCompressibleEulerEquations1D(eos)
    @test Trixi.get_name(equations) ==
          "NonIdealCompressibleEulerEquations1D{VanDerWaals}"
    q = SVector(2.0, 0.1, 10.0)
    V, v1, T = q
    u = thermo2cons(q, equations)

    @test density(u, equations) ≈ 0.5
    @test velocity(u, equations) ≈ 0.1
    @test density_pressure(u, equations) ≈ u[1] * pressure(V, T, eos)
    @test energy_internal_specific(u, equations) ≈ energy_internal_specific(V, T, eos)

    @test ForwardDiff.gradient(u -> entropy(u, equations), u) ≈
          cons2entropy(u, equations)
    @test flux_lax_friedrichs(u, u, 1, equations) ≈ flux(u, 1, equations)
    @test flux_hll(u, u, 1, equations) ≈ flux(u, 1, equations)

    @test flux_terashima_etal(u, u, 1, equations) ≈ flux(u, 1, equations)
    @test flux_central_terashima_etal(u, u, 1, equations) ≈ flux(u, 1, equations)

    # check that the fallback temperature and specialized temperature
    # return the same value
    V, v1, T = cons2thermo(u, equations)
    e_internal = energy_internal_specific(V, T, eos)
    @test temperature(V, e_internal, eos) ≈
          invoke(temperature, Tuple{Any, Any, Trixi.AbstractEquationOfState}, V,
                 e_internal, eos)
    @test cons2prim(u, equations) ≈
          SVector(u[1], v1, pressure(u, equations))

    # check that fallback calc_pressure_derivatives matches specialized routines
    @test Trixi.calc_pressure_derivatives(V, T, eos)[1] ≈
          invoke(Trixi.calc_pressure_derivatives,
                 Tuple{Any, Any, Trixi.AbstractEquationOfState}, V, T, eos)[1]
    @test Trixi.calc_pressure_derivatives(V, T, eos)[2] ≈
          invoke(Trixi.calc_pressure_derivatives,
                 Tuple{Any, Any, Trixi.AbstractEquationOfState}, V, T, eos)[2]

    eos = ThermallyPerfectGas9PolyFit()

    equations = NonIdealCompressibleEulerEquations1D(eos)

    # Mach 26 at 120km altitude, data from US Standard Atmosphere 1976
    rho = 2.22e-8 # [kg/m^3]
    V = inv(rho) # [m^3/kg]

    a = 380.4 # [m/s]
    v1 = 26 * a # [m/s]
    T = 360 # [K]
    q = SVector(V, v1, T)
    u = thermo2cons(q, equations)

    @test flux_lax_friedrichs(u, u, 1, equations) ≈ flux(u, 1, equations)
    @test flux_hll(u, u, 1, equations) ≈ flux(u, 1, equations)

    @test flux_terashima_etal(u, u, 1, equations) ≈ flux(u, 1, equations)
    @test flux_central_terashima_etal(u, u, 1, equations) ≈ flux(u, 1, equations)
end

@testitem "Unit: Test consistency (fluxes, entropy/cons2entropy) for NonIdealCompressibleEulerEquations2D" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    eos = VanDerWaals(; a = 10, b = 0.01, R = 287, gamma = 1.4)
    equations = NonIdealCompressibleEulerEquations2D(eos)
    q = SVector(2.0, 0.1, 0.2, 10.0)
    V, v1, v2, T = q
    u = thermo2cons(q, equations)

    @test density(u, equations) ≈ 0.5
    @test velocity(u, equations) ≈ SVector(0.1, 0.2)
    @test velocity(u, 1, equations) ≈ 0.1
    @test density_pressure(u, equations) ≈ u[1] * pressure(V, T, eos)
    @test energy_internal_specific(u, equations) ≈ energy_internal_specific(V, T, eos)
    @test energy_internal_specific(u, equations) ≈ energy_internal(u, equations) * V
    @test ForwardDiff.gradient(u -> entropy(u, equations), u) ≈
          cons2entropy(u, equations)
    for orientation in (1, 2)
        @test flux_lax_friedrichs(u, u, orientation, equations) ≈
              flux(u, orientation, equations)
        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)

        @test flux_terashima_etal(u, u, orientation, equations) ≈
              flux(u, orientation, equations)
        @test flux_central_terashima_etal(u, u, orientation, equations) ≈
              flux(u, orientation, equations)
    end

    normal_direction = SVector(1, 2) / norm(SVector(1, 2))
    @test flux(u, normal_direction, equations) ≈
          flux(u, 1, equations) * normal_direction[1] +
          flux(u, 2, equations) * normal_direction[2]

    u_ll = u
    u_rr = thermo2cons(SVector(2.5, 0.2, 0.1, 8.0), equations)
    @test flux_terashima_etal(u_ll, u_rr, normal_direction, equations) ≈
          flux_terashima_etal(u_ll, u_rr, 1, equations) * normal_direction[1] +
          flux_terashima_etal(u_ll, u_rr, 2, equations) * normal_direction[2]
    @test flux_central_terashima_etal(u_ll, u_rr, normal_direction, equations) ≈
          flux_central_terashima_etal(u_ll, u_rr, 1, equations) * normal_direction[1] +
          flux_central_terashima_etal(u_ll, u_rr, 2, equations) * normal_direction[2]

    for _flux_function in (flux_lax_friedrichs, min_max_speed_davis)
        @test all(_flux_function(u_ll, u_rr, 1, equations) .≈
                  _flux_function(u_ll, u_rr, SVector(1, 0), equations))
        @test all(_flux_function(u_ll, u_rr, 2, equations) .≈
                  _flux_function(u_ll, u_rr, SVector(0, 1), equations))
    end

    # check consistency of slip wall boundary conditions
    for orientation in (1, 2)
        x, t = 0, 0
        direction = 1 # this variable is not used in `boundary_condition_slip_wall`
        normal_direction = orientation == 1 ? SVector(1.0, 0.0) : SVector(0.0, 1.0)
        @test boundary_condition_slip_wall(u, orientation, direction, x, t,
                                           flux_lax_friedrichs, equations) ≈
              boundary_condition_slip_wall(u, normal_direction, x, t,
                                           flux_lax_friedrichs, equations)
    end

    # check that the fallback temperature and specialized temperature
    # return the same value
    V, v1, v2, T = cons2thermo(u, equations)
    e = energy_internal_specific(V, T, eos)
    @test temperature(V, e, eos) ≈
          invoke(temperature, Tuple{Any, Any, Trixi.AbstractEquationOfState}, V, e, eos)

    # check that fallback calc_pressure_derivatives matches specialized routines
    @test Trixi.calc_pressure_derivatives(V, T, eos)[1] ≈
          invoke(Trixi.calc_pressure_derivatives,
                 Tuple{Any, Any, Trixi.AbstractEquationOfState}, V, T, eos)[1]
    @test Trixi.calc_pressure_derivatives(V, T, eos)[2] ≈
          invoke(Trixi.calc_pressure_derivatives,
                 Tuple{Any, Any, Trixi.AbstractEquationOfState}, V, T, eos)[2]

    eos = ThermallyPerfectGas9PolyFit()

    equations = NonIdealCompressibleEulerEquations2D(eos)

    # Mach 26 at 120km altitude, data from US Standard Atmosphere 1976
    rho = 2.22e-8 # [kg/m^3]
    V = inv(rho) # [m^3/kg]

    a = 380.4 # [m/s]

    aoa = deg2rad(40.0) # angle of attack

    v1 = 26 * a * cos(aoa) # [m/s]
    v2 = 26 * a * sin(aoa) # [m/s]
    T = 360 # [K]
    q = SVector(V, v1, v2, T)
    u = thermo2cons(q, equations)

    for orientation in (1, 2)
        @test flux_lax_friedrichs(u, u, orientation, equations) ≈
              flux(u, orientation, equations)
        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)

        @test flux_terashima_etal(u, u, orientation, equations) ≈
              flux(u, orientation, equations)
        @test flux_central_terashima_etal(u, u, orientation, equations) ≈
              flux(u, orientation, equations)
    end
end

@testitem "Unit: StepsizeCallback" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Ensure a proper error is thrown if used with adaptive time integration schemes
    @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem",
                                 "elixir_advection_diffusion.jl"),
                        tspan=(0, 0.05))

    @test_throws ArgumentError solve(ode, alg; ode_default_options()...,
                                     callback = StepsizeCallback(cfl = 1.0))
end

@testitem "Unit: TimeSeriesCallback" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Test the 2D TreeMesh version of the callback and some warnings
    @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem",
                                 "elixir_acoustics_gaussian_source.jl"),
                        tspan=(0, 0.05))

    point_data_1 = time_series.affect!.point_data[1]
    @test all(isapprox.(point_data_1[1:7],
                        [-2.4417734981719132e-5, -3.4296207289200194e-5,
                            0.0018130846385739788, -0.5, 0.25, 1.0, 1.0]))
    @test_throws DimensionMismatch Trixi.get_elements_by_coordinates!([1, 2],
                                                                      rand(2, 4), mesh,
                                                                      solver, nothing)
    @test_nowarn show(stdout, time_series)
    @test_throws ArgumentError TimeSeriesCallback(semi, [(1.0, 1.0)]; interval = -1)
    @test_throws ArgumentError TimeSeriesCallback(semi, [1.0 1.0 1.0; 2.0 2.0 2.0])
end

@testitem "Unit: resize! RelaxationIntegrators" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    equations = LinearScalarAdvectionEquation1D(42.0)
    solver = DGSEM(polydeg = 0, surface_flux = flux_ranocha)
    mesh = TreeMesh((0.0,), (1.0,),
                    initial_refinement_level = 2,
                    periodicity = true)
    semi = SemidiscretizationHyperbolic(mesh, equations,
                                        initial_condition_convergence_test,
                                        solver;
                                        boundary_conditions = boundary_condition_periodic)
    u0 = zeros(4)
    tspan = (0.0, 1.0)
    ode = semidiscretize(semi, tspan)

    ode_alg = Trixi.RelaxationRK44() # SubDiagonalAlgorithm
    integrator = Trixi.init(ode, ode_alg; dt = 1.0) # SubDiagonalRelaxationIntegrator

    resize!(integrator, 1001)
    @test length(integrator.u) == 1001
    @test length(integrator.du) == 1001
    @test length(integrator.u_tmp) == 1001
    @test length(integrator.direction) == 1001

    ode_alg = Trixi.RelaxationCKL54() # vanderHouwenAlgorithm
    integrator = Trixi.init(ode, ode_alg; dt = 1.0) # vanderHouwenRelaxationIntegrator

    resize!(integrator, 42)
    @test length(integrator.u) == 42
    @test length(integrator.du) == 42
    @test length(integrator.u_tmp) == 42
    @test length(integrator.k_prev) == 42
    @test length(integrator.direction) == 42
end

@testitem "Unit: Consistency check for single point flux: CEMCE" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4),
                                                           gas_constants = (0.4, 0.4))
    u = SVector(0.1, -0.5, 1.0, 1.0, 2.0)

    orientations = [1, 2]
    for orientation in orientations
        @test flux(u, orientation, equations) ≈
              flux_ranocha(u, u, orientation, equations)
    end
end

@testitem "Unit: Consistency check for HLL flux (naive): CEE" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    flux_hll = FluxHLL(min_max_speed_naive)

    # Set up equations and dummy conservative variables state
    equations = CompressibleEulerEquations1D(1.4)
    u = SVector(1.1, 2.34, 5.5)

    orientations = [1]
    for orientation in orientations
        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end

    equations = CompressibleEulerEquations2D(1.4)
    u = SVector(1.1, -0.5, 2.34, 5.5)

    orientations = [1, 2]
    for orientation in orientations
        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end

    equations = CompressibleEulerEquations3D(1.4)
    u = SVector(1.1, -0.5, 2.34, 2.4, 5.5)

    orientations = [1, 2, 3]
    for orientation in orientations
        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end
end

@testitem "Unit: Consistency check for flux_chan_etal: CEEQ" setup=[Setup, UnitTests] tags=[:misc_part1] begin

    # Set up equations and dummy conservative variables state
    equations = CompressibleEulerEquationsQuasi1D(1.4)
    u = SVector(1.1, 2.34, 5.5, 2.73)

    orientations = [1]
    for orientation in orientations
        @test flux_chan_etal(u, u, orientation, equations) ≈
              flux(u, orientation, equations)
    end
end

@testitem "Unit: Consistency check for HLL flux (naive): LEE" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    flux_hll = FluxHLL(min_max_speed_naive)

    equations = LinearizedEulerEquations2D(SVector(1.0, 1.0), 1.0, 1.0)
    u = SVector(1.1, -0.5, 2.34, 5.5)

    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

@testitem "Unit: Consistency check for HLL flux (naive): MHD" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    flux_hll = FluxHLL(min_max_speed_naive)

    equations = IdealGlmMhdEquations1D(1.4)
    u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1),
        SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2)]

    for u in u_values
        @test flux_hll(u, u, 1, equations) ≈ flux(u, 1, equations)
    end

    equations = IdealGlmMhdEquations2D(1.4, 5.0) #= c_h =#
    normal_directions = [SVector(1.0, 0.0),
        SVector(0.0, 1.0),
        SVector(0.5, -0.5),
        SVector(-1.2, 0.3)]
    orientations = [1, 2]

    u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1, 0.0),
        SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2)]

    for u in u_values, orientation in orientations
        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end

    for u in u_values, normal_direction in normal_directions
        @test flux_hll(u, u, normal_direction, equations) ≈
              flux(u, normal_direction, equations)
    end

    equations = IdealGlmMhdEquations3D(1.4, 5.0) #= c_h =#
    normal_directions = [SVector(1.0, 0.0, 0.0),
        SVector(0.0, 1.0, 0.0),
        SVector(0.0, 0.0, 1.0),
        SVector(0.5, -0.5, 0.2),
        SVector(-1.2, 0.3, 1.4)]
    orientations = [1, 2, 3]

    u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1, 0.0),
        SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2)]

    for u in u_values, orientation in orientations
        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end

    for u in u_values, normal_direction in normal_directions
        @test flux_hll(u, u, normal_direction, equations) ≈
              flux(u, normal_direction, equations)
    end
end

@testitem "Unit: Consistency check for HLL flux with Davis wave speed estimates: CEE" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    flux_hll = FluxHLL(min_max_speed_davis)

    # Set up equations and dummy conservative variables state
    equations = CompressibleEulerEquations1D(1.4)
    u = SVector(1.1, 2.34, 5.5)

    orientations = [1]
    for orientation in orientations
        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end

    equations = CompressibleEulerEquations2D(1.4)
    u = SVector(1.1, -0.5, 2.34, 5.5)

    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

    equations = CompressibleEulerEquations3D(1.4)
    u = SVector(1.1, -0.5, 2.34, 2.4, 5.5)

    orientations = [1, 2, 3]
    for orientation in orientations
        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end

    normal_directions = [SVector(1.0, 0.0, 0.0),
        SVector(0.0, 1.0, 0.0),
        SVector(0.0, 0.0, 1.0),
        SVector(0.5, -0.5, 0.2),
        SVector(-1.2, 0.3, 1.4)]

    for normal_direction in normal_directions
        @test flux_hll(u, u, normal_direction, equations) ≈
              flux(u, normal_direction, equations)
    end
end

@testitem "Unit: Consistency check for HLL flux with Davis wave speed estimates: Polytropic CEE" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] 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

@testitem "Unit: Consistency check for Winters flux: Polytropic CEE" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] 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

@testitem "Unit: Consistency check for Lax-Friedrich flux: Polytropic CEE" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] 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_lax_friedrichs(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_lax_friedrichs(u, u, normal_direction, equations) ≈
                  flux(u, normal_direction, equations)
        end
    end
end

@testitem "Unit: Consistency check for HLL flux with Davis wave speed estimates: LEE" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    flux_hll = FluxHLL(min_max_speed_davis)

    equations = LinearizedEulerEquations2D(SVector(1.0, 1.0), 1.0, 1.0)
    u = SVector(1.1, -0.5, 2.34, 5.5)

    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

@testitem "Unit: Consistency check for HLL flux with Davis wave speed estimates: MHD" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    flux_hll = FluxHLL(min_max_speed_davis)

    equations = IdealGlmMhdEquations1D(1.4)
    u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1),
        SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2)]

    for u in u_values
        @test flux_hll(u, u, 1, equations) ≈ flux(u, 1, equations)
    end

    equations = IdealGlmMhdEquations2D(1.4, 5.0) #= c_h =#
    normal_directions = [SVector(1.0, 0.0),
        SVector(0.0, 1.0),
        SVector(0.5, -0.5),
        SVector(-1.2, 0.3)]
    orientations = [1, 2]

    u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1, 0.0),
        SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2)]

    for u in u_values, orientation in orientations
        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end

    for u in u_values, normal_direction in normal_directions
        @test flux_hll(u, u, normal_direction, equations) ≈
              flux(u, normal_direction, equations)
    end

    equations = IdealGlmMhdEquations3D(1.4, 5.0) #= c_h =#
    normal_directions = [SVector(1.0, 0.0, 0.0),
        SVector(0.0, 1.0, 0.0),
        SVector(0.0, 0.0, 1.0),
        SVector(0.5, -0.5, 0.2),
        SVector(-1.2, 0.3, 1.4)]
    orientations = [1, 2, 3]

    u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1, 0.0),
        SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2)]

    for u in u_values, orientation in orientations
        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end

    for u in u_values, normal_direction in normal_directions
        @test flux_hll(u, u, normal_direction, equations) ≈
              flux(u, normal_direction, equations)
    end
end

@testitem "Unit: Consistency check for HLLE flux: CEE" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Set up equations and dummy conservative variables state
    equations = CompressibleEulerEquations1D(1.4)
    u = SVector(1.1, 2.34, 5.5)

    orientations = [1]
    for orientation in orientations
        @test flux_hlle(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end

    equations = CompressibleEulerEquations2D(1.4)
    u = SVector(1.1, -0.5, 2.34, 5.5)

    orientations = [1, 2]
    for orientation in orientations
        @test flux_hlle(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_hlle(u, u, normal_direction, equations) ≈
              flux(u, normal_direction, equations)
    end

    equations = CompressibleEulerEquations3D(1.4)
    u = SVector(1.1, -0.5, 2.34, 2.4, 5.5)

    orientations = [1, 2, 3]
    for orientation in orientations
        @test flux_hlle(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end

    normal_directions = [SVector(1.0, 0.0, 0.0),
        SVector(0.0, 1.0, 0.0),
        SVector(0.0, 0.0, 1.0),
        SVector(0.5, -0.5, 0.2),
        SVector(-1.2, 0.3, 1.4)]

    for normal_direction in normal_directions
        @test flux_hlle(u, u, normal_direction, equations) ≈
              flux(u, normal_direction, equations)
    end
end

@testitem "Unit: Consistency check for HLLE flux: MHD" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    equations = IdealGlmMhdEquations1D(1.4)
    u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1),
        SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2)]

    for u in u_values
        @test flux_hlle(u, u, 1, equations) ≈ flux(u, 1, equations)
        @test flux_hllc(u, u, 1, equations) ≈ flux(u, 1, equations)
    end

    equations = IdealGlmMhdEquations2D(1.4, 5.0) #= c_h =#
    normal_directions = [SVector(1.0, 0.0),
        SVector(0.0, 1.0),
        SVector(0.5, -0.5),
        SVector(-1.2, 0.3)]
    orientations = [1, 2]

    u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1, 0.0),
        SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2)]

    for u in u_values, orientation in orientations
        @test flux_hlle(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end

    for u in u_values, normal_direction in normal_directions
        @test flux_hlle(u, u, normal_direction, equations) ≈
              flux(u, normal_direction, equations)
    end

    equations = IdealGlmMhdEquations3D(1.4, 5.0) #= c_h =#
    normal_directions = [SVector(1.0, 0.0, 0.0),
        SVector(0.0, 1.0, 0.0),
        SVector(0.0, 0.0, 1.0),
        SVector(0.5, -0.5, 0.2),
        SVector(-1.2, 0.3, 1.4)]
    orientations = [1, 2, 3]

    u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1, 0.0),
        SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2)]

    for u in u_values, orientation in orientations
        @test flux_hlle(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end

    for u in u_values, normal_direction in normal_directions
        @test flux_hlle(u, u, normal_direction, equations) ≈
              flux(u, normal_direction, equations)
    end
end

@testitem "Unit: Consistency check for HLLC flux: CEE" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Set up equations and dummy conservative variables state
    equations = CompressibleEulerEquations2D(1.4)
    u = SVector(1.1, -0.5, 2.34, 5.5)

    orientations = [1, 2]
    for orientation in orientations
        @test flux_hllc(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_hllc(u, u, normal_direction, equations) ≈
              flux(u, normal_direction, equations)
    end

    # check consistency between 1D and 2D HLLC fluxes
    u_1d = SVector(1.1, -0.5, 5.5)
    u_2d = SVector(u_1d[1], u_1d[2], 0.0, u_1d[3])
    normal_1d = SVector(-0.3)
    normal_2d = SVector(normal_1d[1], 0.0)
    equations_1d = CompressibleEulerEquations1D(1.4)
    equations_2d = CompressibleEulerEquations2D(1.4)
    flux_1d = flux_hllc(u_1d, u_1d, normal_1d, equations_1d)
    flux_2d = flux_hllc(u_2d, u_2d, normal_2d, equations_2d)
    @test flux_1d ≈ flux(u_1d, normal_1d, equations_1d)
    @test flux_1d ≈ flux_2d[[1, 2, 4]]

    # test when u_ll is not the same as u_rr
    u_rr_1d = SVector(2.1, 0.3, 0.1)
    u_rr_2d = SVector(u_rr_1d[1], u_rr_1d[2], 0.0, u_rr_1d[3])
    flux_1d = flux_hllc(u_1d, u_rr_1d, normal_1d, equations_1d)
    flux_2d = flux_hllc(u_2d, u_rr_2d, normal_2d, equations_2d)
    @test flux_1d ≈ flux_2d[[1, 2, 4]]

    equations = CompressibleEulerEquations3D(1.4)
    u = SVector(1.1, -0.5, 2.34, 2.4, 5.5)

    orientations = [1, 2, 3]
    for orientation in orientations
        @test flux_hllc(u, u, orientation, equations) ≈ flux(u, orientation, equations)
    end

    normal_directions = [SVector(1.0, 0.0, 0.0),
        SVector(0.0, 1.0, 0.0),
        SVector(0.0, 0.0, 1.0),
        SVector(0.5, -0.5, 0.2),
        SVector(-1.2, 0.3, 1.4)]

    for normal_direction in normal_directions
        @test flux_hllc(u, u, normal_direction, equations) ≈
              flux(u, normal_direction, equations)
    end
end

@testitem "Unit: Consistency check for Godunov flux" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Set up equations and dummy conservative variables state
    # Burgers' Equation

    equation = InviscidBurgersEquation1D()
    u_values = [SVector(42.0), SVector(-42.0)]

    orientations = [1]
    for orientation in orientations, u in u_values
        @test flux_godunov(u, u, orientation, equation) ≈ flux(u, orientation, equation)
    end

    # Linear Advection 1D
    equation = LinearScalarAdvectionEquation1D(-4.2)
    u = SVector(3.14159)

    orientations = [1]
    for orientation in orientations
        @test flux_godunov(u, u, orientation, equation) ≈ flux(u, orientation, equation)
    end

    # Linear Advection 2D
    equation = LinearScalarAdvectionEquation2D(-4.2, 2.4)
    u = SVector(3.14159)

    orientations = [1, 2]
    for orientation in orientations
        @test flux_godunov(u, u, orientation, equation) ≈ flux(u, orientation, equation)
    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_godunov(u, u, normal_direction, equation) ≈
              flux(u, normal_direction, equation)
    end

    # Linear Advection 3D
    equation = LinearScalarAdvectionEquation3D(-4.2, 2.4, 1.2)
    u = SVector(3.14159)

    orientations = [1, 2, 3]
    for orientation in orientations
        @test flux_godunov(u, u, orientation, equation) ≈ flux(u, orientation, equation)
    end

    normal_directions = [SVector(1.0, 0.0, 0.0),
        SVector(0.0, 1.0, 0.0),
        SVector(0.0, 0.0, 1.0),
        SVector(0.5, -0.5, 0.2),
        SVector(-1.2, 0.3, 1.4)]

    for normal_direction in normal_directions
        @test flux_godunov(u, u, normal_direction, equation) ≈
              flux(u, normal_direction, equation)
    end

    # Linearized Euler 2D
    equation = LinearizedEulerEquations2D(v_mean_global = (0.5, -0.7),
                                          c_mean_global = 1.1,
                                          rho_mean_global = 1.2)
    u_values = [SVector(1.0, 0.5, -0.7, 1.0),
        SVector(1.5, -0.2, 0.1, 5.0)]

    orientations = [1, 2]
    for orientation in orientations, u in u_values
        @test flux_godunov(u, u, orientation, equation) ≈ flux(u, orientation, equation)
    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, u in u_values
        @test flux_godunov(u, u, normal_direction, equation) ≈
              flux(u, normal_direction, equation)
    end
end

@testitem "Unit: Consistency check for entropy-conserving Burgers flux" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    equations = InviscidBurgersEquation1D()
    u_ll = SVector(2.0)
    u_rr = SVector(-1.0)

    for normal_direction in (SVector(1.0), SVector(-1.2))
        @test flux_ec(u_ll, u_rr, normal_direction, equations) ≈
              normal_direction[1] * flux_ec(u_ll, u_rr, 1, equations)
    end

    u = SVector(42.0)
    normal_direction = SVector(-1.2)
    @test flux_ec(u, u, normal_direction, equations) ≈
          flux(u, normal_direction, equations)
end

@testitem "Unit: Consistency check for Engquist-Osher flux" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Set up equations and dummy conservative variables state
    equation = InviscidBurgersEquation1D()
    u_values = [SVector(42.0), SVector(-42.0)]

    orientations = [1]
    for orientation in orientations, u in u_values
        @test Trixi.flux_engquist_osher(u, u, orientation, equation) ≈
              flux(u, orientation, equation)
    end

    equation = LinearScalarAdvectionEquation1D(-4.2)
    u = SVector(3.14159)

    orientations = [1]
    for orientation in orientations
        @test Trixi.flux_engquist_osher(u, u, orientation, equation) ≈
              flux(u, orientation, equation)
    end
end

@testitem "Unit: Flux consistency checks LinearElasticityEquations1D" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    rho = 7800.0 # kg/m³
    lambda = 9.3288e10
    mu = lambda
    equations = LinearElasticityEquations1D(rho = rho, mu = mu, lambda = lambda)

    u = SVector(1.42, 2.666)

    orientation = 1
    @test flux_central(u, u, orientation, equations) ≈
          flux(u, orientation, equations)

    @test flux_lax_friedrichs(u, u, orientation, equations) ≈
          flux(u, orientation, equations)

    @test flux_hll(u, u, orientation, equations) ≈
          flux(u, orientation, equations)
end

@testitem "Unit: Consistency check for `gradient_conservative` routine" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    # Set up conservative variables, equations
    u = [
        0.5011914484393387,
        0.8829127712445113,
        0.43024132987932817,
        0.7560616633050348
    ]

    equations = CompressibleEulerEquations2D(1.4)

    # Define wrapper function for pressure in order to call default implementation
    function pressure_test(u, equations)
        return pressure(u, equations)
    end

    @test Trixi.gradient_conservative(pressure_test, u, equations) ≈
          Trixi.gradient_conservative(pressure, u, equations)
end

@testitem "Unit: Equivalent Fluxes" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # Set up equations and dummy conservative variables state
    # Burgers' Equation

    equation = InviscidBurgersEquation1D()
    u_values = [SVector(42.0), SVector(-42.0)]

    orientations = [1]
    for orientation in orientations, u in u_values
        @test flux_godunov(0.75 * u, u, orientation, equation) ≈
              Trixi.flux_engquist_osher(0.75 * u, u, orientation, equation)
    end

    # Linear Advection 1D
    equation = LinearScalarAdvectionEquation1D(-4.2)
    u = SVector(3.14159)

    orientations = [1]
    for orientation in orientations
        @test flux_godunov(0.5 * u, u, orientation, equation) ≈
              flux_lax_friedrichs(0.5 * u, u, orientation, equation)
        @test flux_godunov(2 * u, u, orientation, equation) ≈
              Trixi.flux_engquist_osher(2 * u, u, orientation, equation)
    end

    # Linear Advection 2D
    equation = LinearScalarAdvectionEquation2D(-4.2, 2.4)
    u = SVector(3.14159)

    orientations = [1, 2]
    for orientation in orientations
        @test flux_godunov(0.25 * u, u, orientation, equation) ≈
              flux_lax_friedrichs(0.25 * u, u, orientation, equation)
    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_godunov(3 * u, u, normal_direction, equation) ≈
              flux_lax_friedrichs(3 * u, u, normal_direction, equation)
    end

    # Linear Advection 3D
    equation = LinearScalarAdvectionEquation3D(-4.2, 2.4, 1.2)
    u = SVector(3.14159)

    orientations = [1, 2, 3]
    for orientation in orientations
        @test flux_godunov(1.5 * u, u, orientation, equation) ≈
              flux_lax_friedrichs(1.5 * u, u, orientation, equation)
    end

    normal_directions = [SVector(1.0, 0.0, 0.0),
        SVector(0.0, 1.0, 0.0),
        SVector(0.0, 0.0, 1.0),
        SVector(0.5, -0.5, 0.2),
        SVector(-1.2, 0.3, 1.4)]

    for normal_direction in normal_directions
        @test flux_godunov(1.3 * u, u, normal_direction, equation) ≈
              flux_lax_friedrichs(1.3 * u, u, normal_direction, equation)
    end
end

@testitem "Unit: Consistency check for LMARS flux" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    equations = CompressibleEulerEquations2D(1.4)
    flux_lmars = FluxLMARS(340)

    normal_directions = [SVector(1.0, 0.0),
        SVector(0.0, 1.0),
        SVector(0.5, -0.5),
        SVector(-1.2, 0.3)]
    orientations = [1, 2]
    u_values = [SVector(1.0, 0.5, -0.7, 1.0),
        SVector(1.5, -0.2, 0.1, 5.0)]

    for u in u_values, orientation in orientations
        @test flux_lmars(u, u, orientation, equations) ≈
              flux(u, orientation, equations)
    end

    for u in u_values, normal_direction in normal_directions
        @test flux_lmars(u, u, normal_direction, equations) ≈
              flux(u, normal_direction, equations)
    end

    equations = CompressibleEulerEquations3D(1.4)
    normal_directions = [SVector(1.0, 0.0, 0.0),
        SVector(0.0, 1.0, 0.0),
        SVector(0.0, 0.0, 1.0),
        SVector(0.5, -0.5, 0.2),
        SVector(-1.2, 0.3, 1.4)]
    orientations = [1, 2, 3]
    u_values = [SVector(1.0, 0.5, -0.7, 0.1, 1.0),
        SVector(1.5, -0.2, 0.1, 0.2, 5.0)]

    for u in u_values, orientation in orientations
        @test flux_lmars(u, u, orientation, equations) ≈
              flux(u, orientation, equations)
    end

    for u in u_values, normal_direction in normal_directions
        @test flux_lmars(u, u, normal_direction, equations) ≈
              flux(u, normal_direction, equations)
    end
end

@testitem "Unit: FluxRotated vs. direct implementation" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    @timed_testset "CompressibleEulerMulticomponentEquations2D" begin
        equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4),
                                                               gas_constants = (0.4,
                                                                                0.4))
        normal_directions = [SVector(1.0, 0.0),
            SVector(0.0, 1.0),
            SVector(0.5, -0.5),
            SVector(-1.2, 0.3)]
        u_values = [SVector(0.1, -0.5, 1.0, 1.0, 2.0),
            SVector(-0.1, -0.3, 1.2, 1.3, 1.4)]

        f_std = flux
        f_rot = FluxRotated(f_std)
        println(typeof(f_std))
        println(typeof(f_rot))
        for u in u_values,
            normal_direction in normal_directions

            @test f_rot(u, normal_direction, equations) ≈
                  f_std(u, normal_direction, equations)
        end
    end

    @timed_testset "CompressibleEulerEquations2D" begin
        equations = CompressibleEulerEquations2D(1.4)
        normal_directions = [SVector(1.0, 0.0),
            SVector(0.0, 1.0),
            SVector(0.5, -0.5),
            SVector(-1.2, 0.3)]
        u_values = [SVector(1.0, 0.5, -0.7, 1.0),
            SVector(1.5, -0.2, 0.1, 5.0)]
        fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber,
            FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle,
            flux_hllc, flux_chandrashekar
        ]

        for f_std in fluxes
            f_rot = FluxRotated(f_std)
            for u_ll in u_values, u_rr in u_values,
                normal_direction in normal_directions

                @test f_rot(u_ll, u_rr, normal_direction, equations) ≈
                      f_std(u_ll, u_rr, normal_direction, equations)
            end
        end
    end

    @timed_testset "CompressibleEulerEquations3D" begin
        equations = CompressibleEulerEquations3D(1.4)
        normal_directions = [SVector(1.0, 0.0, 0.0),
            SVector(0.0, 1.0, 0.0),
            SVector(0.0, 0.0, 1.0),
            SVector(0.5, -0.5, 0.2),
            SVector(-1.2, 0.3, 1.4)]
        u_values = [SVector(1.0, 0.5, -0.7, 0.1, 1.0),
            SVector(1.5, -0.2, 0.1, 0.2, 5.0)]
        fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber,
            FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle,
            flux_hllc, flux_chandrashekar
        ]

        for f_std in fluxes
            f_rot = FluxRotated(f_std)
            for u_ll in u_values, u_rr in u_values,
                normal_direction in normal_directions

                @test f_rot(u_ll, u_rr, normal_direction, equations) ≈
                      f_std(u_ll, u_rr, normal_direction, equations)
            end
        end
    end

    @timed_testset "IdealGlmMhdEquations2D" begin
        equations = IdealGlmMhdEquations2D(1.4, 5.0) #= c_h =#
        normal_directions = [SVector(1.0, 0.0),
            SVector(0.0, 1.0),
            SVector(0.5, -0.5),
            SVector(-1.2, 0.3)]
        u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1, 0.0),
            SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2)]
        fluxes = [
            flux_central,
            flux_hindenlang_gassner,
            FluxHLL(min_max_speed_davis),
            flux_hlle
        ]

        for f_std in fluxes
            f_rot = FluxRotated(f_std)
            for u_ll in u_values, u_rr in u_values,
                normal_direction in normal_directions

                @test f_rot(u_ll, u_rr, normal_direction, equations) ≈
                      f_std(u_ll, u_rr, normal_direction, equations)
            end
        end
    end

    @timed_testset "IdealGlmMhdEquations3D" begin
        equations = IdealGlmMhdEquations3D(1.4, 5.0) #= c_h =#
        normal_directions = [SVector(1.0, 0.0, 0.0),
            SVector(0.0, 1.0, 0.0),
            SVector(0.0, 0.0, 1.0),
            SVector(0.5, -0.5, 0.2),
            SVector(-1.2, 0.3, 1.4)]
        u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1, 0.0),
            SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2)]
        fluxes = [
            flux_central,
            flux_hindenlang_gassner,
            FluxHLL(min_max_speed_davis),
            flux_hlle
        ]

        for f_std in fluxes
            f_rot = FluxRotated(f_std)
            for u_ll in u_values, u_rr in u_values,
                normal_direction in normal_directions

                @test f_rot(u_ll, u_rr, normal_direction, equations) ≈
                      f_std(u_ll, u_rr, normal_direction, equations)
            end
        end
    end
end

@testitem "Unit: DissipationMatrixWintersEtal entropy dissipation and consistency tests" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    equations = CompressibleEulerEquations1D(1.4)
    dissipation_matrix_winters_etal = DissipationMatrixWintersEtal()

    # test constant preservation and entropy dissipation vector
    u_ll = prim2cons(SVector(1, 0, 2.0), equations)
    u_rr = prim2cons(SVector(1.1, 0, 2.0), equations)
    v_ll = cons2entropy(u_ll, equations)
    v_rr = cons2entropy(u_rr, equations)
    @test norm(dissipation_matrix_winters_etal(u_ll, u_rr, SVector(1.0), equations)) <
          100 * eps()
    @test dot(v_ll - v_rr,
              dissipation_matrix_winters_etal(u_ll, u_rr, SVector(1.0), equations)) ≥ 0

    # test non-unit vector
    u_ll = prim2cons(SVector(rand(), randn(), rand()), equations)
    u_rr = prim2cons(SVector(rand(), randn(), rand()), equations)
    v_ll = cons2entropy(u_ll, equations)
    v_rr = cons2entropy(u_rr, equations)
    @test dot(v_ll - v_rr,
              dissipation_matrix_winters_etal(u_ll, u_rr, SVector(1.0), equations)) ≥ 0
    @test dissipation_matrix_winters_etal(u_ll, u_rr, SVector(0.1), equations) ≈
          0.1 * dissipation_matrix_winters_etal(u_ll, u_rr, SVector(1.0), equations)

    equations = CompressibleEulerEquations2D(1.4)

    # test that 2D flux is consistent with 1D matrix flux
    u_ll = prim2cons(SVector(1, 0.1, 0, 1.0), equations)
    u_rr = prim2cons(SVector(1.1, -0.2, 0, 2.0), equations)
    v_ll = cons2entropy(u_ll, equations)
    v_rr = cons2entropy(u_rr, equations)
    normal = SVector(1.0, 0.0)
    ids = [1, 2, 4] # indices of 1D variables/fluxes within the 2D solution
    @test dissipation_matrix_winters_etal(u_ll, u_rr, normal, equations)[ids] ≈
          dissipation_matrix_winters_etal(u_ll[ids], u_rr[ids],
                                          SVector(1.0),
                                          CompressibleEulerEquations1D(1.4))

    # test 2D entropy dissipation
    u_ll = prim2cons(SVector(1, 1, -3, 100.0), equations)
    u_rr = prim2cons(SVector(100, -2, 4, 1.0), equations)
    v_ll = cons2entropy(u_ll, equations)
    v_rr = cons2entropy(u_rr, equations)
    dissipation = dissipation_matrix_winters_etal(u_ll, u_rr, SVector(1.0, 1.0),
                                                  equations)
    @test dot(v_ll - v_rr, dissipation) ≥ 0

    # test non-unit vector
    normal_direction = SVector(1.0, 2.0)
    @test dissipation_matrix_winters_etal(u_ll, u_rr, normal_direction, equations) ≈
          norm(normal_direction) * dissipation_matrix_winters_etal(u_ll, u_rr,
                                          normal_direction / norm(normal_direction),
                                          equations)

    # test that 3D flux is consistent with 1D and 2D versions
    equations = CompressibleEulerEquations3D(1.4)
    dissipation_matrix_winters_etal = DissipationMatrixWintersEtal()

    # test for consistency with 1D and 2D flux
    u_ll = prim2cons(SVector(1, 0.1, 0, 0, 1.0), equations)
    u_rr = prim2cons(SVector(1.1, -0.2, 0, 0, 2.0), equations)
    v_ll = cons2entropy(u_ll, equations)
    v_rr = cons2entropy(u_rr, equations)
    dissipation_3d = dissipation_matrix_winters_etal(u_ll, u_rr, SVector(1.0, 0.0, 0.0),
                                                     equations)
    dissipation_1d = dissipation_matrix_winters_etal(u_ll[[1, 2, 5]], u_rr[[1, 2, 5]],
                                                     SVector(1.0),
                                                     CompressibleEulerEquations1D(1.4))
    @test dissipation_3d[[1, 2, 5]] ≈ dissipation_1d

    u_ll = prim2cons(SVector(1, 0.1, 0.2, 0, 1.0), equations)
    u_rr = prim2cons(SVector(1.1, -0.2, -0.3, 0, 2.0), equations)
    v_ll = cons2entropy(u_ll, equations)
    v_rr = cons2entropy(u_rr, equations)
    dissipation_3d = dissipation_matrix_winters_etal(u_ll, u_rr, SVector(1.0, 1.0, 0.0),
                                                     equations)
    dissipation_2d = dissipation_matrix_winters_etal(u_ll[[1, 2, 3, 5]],
                                                     u_rr[[1, 2, 3, 5]],
                                                     SVector(1.0, 1.0),
                                                     CompressibleEulerEquations2D(1.4))
    @test dissipation_3d[[1, 2, 3, 5]] ≈ dissipation_2d

    # test 3D entropy dissipation
    u_ll = prim2cons(SVector(1, 0.1, 0.2, 0.3, 1.0), equations)
    u_rr = prim2cons(SVector(1.1, -0.2, -0.3, 0.4, 2.0), equations)
    v_ll = cons2entropy(u_ll, equations)
    v_rr = cons2entropy(u_rr, equations)
    dissipation_3d = dissipation_matrix_winters_etal(u_ll, u_rr, SVector(1.0, 2.0, 3.0),
                                                     equations)
    @test dot(v_ll - v_rr, dissipation_3d) ≥ 0

    # test non-unit vector
    normal_direction = SVector(1.0, 2.0, 3.0)
    @test dissipation_matrix_winters_etal(u_ll, u_rr, normal_direction, equations) ≈
          norm(normal_direction) * dissipation_matrix_winters_etal(u_ll, u_rr,
                                          normal_direction / norm(normal_direction),
                                          equations)
end

@testitem "Unit: Equivalent Wave Speed Estimates" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    @timed_testset "Linearized Euler 3D" begin
        equations = LinearizedEulerEquations3D(v_mean_global = (0.42, 0.37, 0.7),
                                               c_mean_global = 1.0,
                                               rho_mean_global = 1.0)

        normal_x = SVector(1.0, 0.0, 0.0)
        normal_y = SVector(0.0, 1.0, 0.0)
        normal_z = SVector(0.0, 0.0, 1.0)

        u_ll = SVector(0.3, 0.5, -0.7, 0.1, 1.0)
        u_rr = SVector(0.5, -0.2, 0.1, 0.2, 5.0)

        @test all(isapprox(x, y)
                  for (x, y) in zip(max_abs_speed_naive(u_ll, u_rr, 1, equations),
                                    max_abs_speed_naive(u_ll, u_rr, normal_x,
                                                        equations)))
        @test all(isapprox(x, y)
                  for (x, y) in zip(max_abs_speed_naive(u_ll, u_rr, 2, equations),
                                    max_abs_speed_naive(u_ll, u_rr, normal_y,
                                                        equations)))
        @test all(isapprox(x, y)
                  for (x, y) in zip(max_abs_speed_naive(u_ll, u_rr, 3, equations),
                                    max_abs_speed_naive(u_ll, u_rr, normal_z,
                                                        equations)))

        @test all(isapprox(x, y)
                  for (x, y) in zip(min_max_speed_naive(u_ll, u_rr, 1, equations),
                                    min_max_speed_naive(u_ll, u_rr, normal_x,
                                                        equations)))
        @test all(isapprox(x, y)
                  for (x, y) in zip(min_max_speed_naive(u_ll, u_rr, 2, equations),
                                    min_max_speed_naive(u_ll, u_rr, normal_y,
                                                        equations)))
        @test all(isapprox(x, y)
                  for (x, y) in zip(min_max_speed_naive(u_ll, u_rr, 3, equations),
                                    min_max_speed_naive(u_ll, u_rr, normal_z,
                                                        equations)))

        @test all(isapprox(x, y)
                  for (x, y) in zip(min_max_speed_davis(u_ll, u_rr, 1, equations),
                                    min_max_speed_davis(u_ll, u_rr, normal_x,
                                                        equations)))
        @test all(isapprox(x, y)
                  for (x, y) in zip(min_max_speed_davis(u_ll, u_rr, 2, equations),
                                    min_max_speed_davis(u_ll, u_rr, normal_y,
                                                        equations)))
        @test all(isapprox(x, y)
                  for (x, y) in zip(min_max_speed_davis(u_ll, u_rr, 3, equations),
                                    min_max_speed_davis(u_ll, u_rr, normal_z,
                                                        equations)))
    end

    @timed_testset "Maxwell 1D" begin
        equations = MaxwellEquations1D()

        u_values_left = [SVector(1.0, 0.0),
            SVector(0.0, 1.0),
            SVector(0.5, -0.5),
            SVector(-1.2, 0.3)]

        u_values_right = [SVector(1.0, 0.0),
            SVector(0.0, 1.0),
            SVector(0.5, -0.5),
            SVector(-1.2, 0.3)]
        for u_ll in u_values_left, u_rr in u_values_right
            @test all(isapprox(x, y)
                      for (x, y) in zip(min_max_speed_naive(u_ll, u_rr, 1, equations),
                                        min_max_speed_davis(u_ll, u_rr, 1, equations)))
        end
    end
end

@testitem "Unit: Equivalent Wave Speed Estimates: max_abs_speed(naive)" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    @timed_testset "AcousticPerturbationEquations2D" begin
        equations = AcousticPerturbationEquations2D(v_mean_global = (0.5, 0.3),
                                                    c_mean_global = 2.0,
                                                    rho_mean_global = 0.9)

        v1_prime_ll_rr = SVector(0.1, 0.2)
        v2_prime_ll_rr = SVector(0.3, 0.4)
        p_prime_scaled_ll_rr = SVector(0.5, 0.6)
        v1_mean_ll_rr = SVector(-0.2, -0.1)
        v2_mean_ll_rr = SVector(-0.9, -1.2)
        c_mean = 2.0 # Same for both to get same wave speed estimates
        rho_mean_ll_rr = SVector(1.3, 1.4)

        u_ll = SVector(v1_prime_ll_rr[1], v2_prime_ll_rr[1], p_prime_scaled_ll_rr[1],
                       v1_mean_ll_rr[1], v2_mean_ll_rr[1], c_mean, rho_mean_ll_rr[1])

        u_rr = SVector(v1_prime_ll_rr[2], v2_prime_ll_rr[2], p_prime_scaled_ll_rr[2],
                       v1_mean_ll_rr[2], v2_mean_ll_rr[2], c_mean, rho_mean_ll_rr[2])

        for orientation in [1, 2]
            @test max_abs_speed_naive(u_ll, u_rr, orientation, equations) ≈
                  max_abs_speed(u_ll, u_rr, orientation, equations)
        end

        normal_directions = [SVector(1.0, 0.0, 0.0),
            SVector(0.0, 1.0, 0.0),
            SVector(0.0, 0.0, 1.0),
            SVector(0.5, -0.5, 0.2),
            SVector(-1.2, 0.3, 1.4)]

        for normal_direction in normal_directions
            @test max_abs_speed_naive(u_ll, u_rr, normal_direction, equations) ≈
                  max_abs_speed(u_ll, u_rr, normal_direction, equations)
        end
    end

    @timed_testset "CompressibleEulerEquations1D" begin
        for gamma in [1.4, 5 / 3, 7 / 5]
            equations = CompressibleEulerEquations1D(gamma)

            p_rho_ratio = 42.0

            rho_ll_rr = SVector(2.0, 1.0)
            v_ll_rr = SVector(0.1, 0.2)
            p_ll_rr = SVector(p_rho_ratio * rho_ll_rr[1], p_rho_ratio * rho_ll_rr[2])

            u_ll = prim2cons(SVector(rho_ll_rr[1], v_ll_rr[1], p_ll_rr[1]), equations)
            u_rr = prim2cons(SVector(rho_ll_rr[2], v_ll_rr[2], p_ll_rr[2]), equations)

            @test max_abs_speed_naive(u_ll, u_rr, 1, equations) ≈
                  max_abs_speed(u_ll, u_rr, 1, equations)
        end
    end

    @timed_testset "Passive tracer equations" begin
        for gamma in [1.4, 5 / 3, 7 / 5]
            flow_equations = CompressibleEulerEquations1D(gamma)
            equations = PassiveTracerEquations(flow_equations, n_tracers = 2)

            p_rho_ratio = 42.0
            xi1_ll, xi1_rr = 0.1, 0.2
            xi2_ll, xi2_rr = 0.3, 0.4

            rho_ll_rr = SVector(2.0, 1.0)
            v_ll_rr = SVector(0.1, 0.2)
            p_ll_rr = SVector(p_rho_ratio * rho_ll_rr[1], p_rho_ratio * rho_ll_rr[2])

            u_ll = prim2cons(SVector(rho_ll_rr[1], v_ll_rr[1], p_ll_rr[1], xi1_ll,
                                     xi2_ll), equations)
            u_rr = prim2cons(SVector(rho_ll_rr[2], v_ll_rr[2], p_ll_rr[2], xi1_rr,
                                     xi2_rr), equations)

            @test max_abs_speed_naive(u_ll, u_rr, 1, equations) ≈
                  max_abs_speed_naive(u_ll, u_rr, 1, flow_equations)
        end
    end

    @timed_testset "CompressibleEulerEquations2D" begin
        for gamma in [1.4, 5 / 3, 7 / 5]
            equations = CompressibleEulerEquations2D(gamma)

            p_rho_ratio = 27.0

            rho_ll_rr = SVector(2.0, 1.0)
            v1_ll_rr = SVector(0.1, 0.2)
            v2_ll_rr = SVector(0.4, 0.3)
            p_ll_rr = SVector(p_rho_ratio * rho_ll_rr[1], p_rho_ratio * rho_ll_rr[2])

            u_ll = prim2cons(SVector(rho_ll_rr[1], v1_ll_rr[1], v2_ll_rr[1],
                                     p_ll_rr[1]), equations)
            u_rr = prim2cons(SVector(rho_ll_rr[2], v1_ll_rr[2], v2_ll_rr[2],
                                     p_ll_rr[2]), equations)

            for orientation in [1, 2]
                @test max_abs_speed_naive(u_ll, u_rr, orientation, equations) ≈
                      max_abs_speed(u_ll, u_rr, 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 max_abs_speed_naive(u_ll, u_rr, normal_direction, equations) ≈
                      max_abs_speed(u_ll, u_rr, normal_direction, equations)
            end
        end
    end

    @timed_testset "CompressibleEulerEquations3D" begin
        for gamma in [1.4, 5 / 3, 7 / 5]
            equations = CompressibleEulerEquations3D(gamma)

            p_rho_ratio = 11.0

            rho_ll_rr = SVector(1.0, 2.0)
            v1_ll_rr = SVector(0.1, 0.2)
            v2_ll_rr = SVector(0.4, 0.3)
            v3_ll_rr = SVector(0.9, 0.8)
            p_ll_rr = SVector(p_rho_ratio * rho_ll_rr[1], p_rho_ratio * rho_ll_rr[2])

            u_ll = prim2cons(SVector(rho_ll_rr[1],
                                     v1_ll_rr[1], v2_ll_rr[1], v3_ll_rr[1],
                                     p_ll_rr[1]), equations)
            u_rr = prim2cons(SVector(rho_ll_rr[2],
                                     v1_ll_rr[2], v2_ll_rr[2], v3_ll_rr[2],
                                     p_ll_rr[2]), equations)

            for orientation in [1, 2, 3]
                @test max_abs_speed_naive(u_ll, u_rr, orientation, equations) ≈
                      max_abs_speed(u_ll, u_rr, orientation, equations)
            end

            normal_directions = [SVector(1.0, 0.0, 0.0),
                SVector(0.0, 1.0, 0.0),
                SVector(0.0, 0.0, 1.0),
                SVector(0.5, -0.5, 0.2),
                SVector(-1.2, 0.3, 1.4)]

            for normal_direction in normal_directions
                @test max_abs_speed_naive(u_ll, u_rr, normal_direction, equations) ≈
                      max_abs_speed(u_ll, u_rr, normal_direction, equations)
            end
        end
    end

    @timed_testset "CompressibleEulerMulticomponentEquations1D" begin
        for gamma in [1.4, 5 / 3, 7 / 5]
            equations = CompressibleEulerMulticomponentEquations1D(gammas = (gamma,
                                                                             gamma),
                                                                   gas_constants = (0.5,
                                                                                    0.4))

            p_rho_ratio = 42.0

            rho1_ll_rr = SVector(2.0, 1.0)
            rho2_ll_rr = SVector(2.0, 1.0)
            v_ll_rr = SVector(0.1, 0.2)
            p_ll_rr = SVector(p_rho_ratio * rho1_ll_rr[1], p_rho_ratio * rho1_ll_rr[2])

            u_ll = prim2cons(SVector(v_ll_rr[1], p_ll_rr[1], rho1_ll_rr[1],
                                     rho2_ll_rr[1]), equations)
            u_rr = prim2cons(SVector(v_ll_rr[2], p_ll_rr[2], rho1_ll_rr[2],
                                     rho2_ll_rr[2]), equations)

            @test max_abs_speed_naive(u_ll, u_rr, 1, equations) ≈
                  max_abs_speed(u_ll, u_rr, 1, equations)
        end
    end

    @timed_testset "CompressibleEulerMulticomponentEquations2D" begin
        for gamma in [1.4, 5 / 3, 7 / 5]
            equations = CompressibleEulerMulticomponentEquations2D(gammas = (gamma,
                                                                             gamma),
                                                                   gas_constants = (0.5,
                                                                                    0.6))

            p_rho_ratio = 27.0

            rho1_ll_rr = SVector(2.0, 1.0)
            rho2_ll_rr = SVector(2.0, 1.0)
            v1_ll_rr = SVector(0.1, 0.2)
            v2_ll_rr = SVector(0.4, 0.3)
            p_ll_rr = SVector(p_rho_ratio * rho1_ll_rr[1], p_rho_ratio * rho1_ll_rr[2])

            u_ll = prim2cons(SVector(v1_ll_rr[1], v2_ll_rr[1], p_ll_rr[1],
                                     rho1_ll_rr[1], rho2_ll_rr[1]), equations)
            u_rr = prim2cons(SVector(v1_ll_rr[2], v2_ll_rr[2], p_ll_rr[2],
                                     rho1_ll_rr[2], rho2_ll_rr[2]), equations)

            for orientation in [1, 2]
                @test max_abs_speed_naive(u_ll, u_rr, orientation, equations) ≈
                      max_abs_speed(u_ll, u_rr, orientation, equations)
            end

            @test max_abs_speed_naive(u_ll, u_rr, 1, equations) ≈
                  max_abs_speed_naive(u_ll, u_rr, SVector(1.0, 0.0), equations)
            @test max_abs_speed_naive(u_ll, u_rr, 2, equations) ≈
                  max_abs_speed_naive(u_ll, u_rr, SVector(0.0, 1.0), equations)
        end
    end

    @timed_testset "CompressibleEulerEquationsQuasi1D" begin
        for gamma in [1.4, 5 / 3, 7 / 5]
            equations = CompressibleEulerEquationsQuasi1D(gamma)

            p_rho_ratio = 11.0

            rho_ll_rr = SVector(1.0, 2.0)
            v_ll_rr = SVector(0.1, 0.2)
            a_ll_rr = SVector(0.3, 0.4)
            p_ll_rr = SVector(p_rho_ratio * rho_ll_rr[1], p_rho_ratio * rho_ll_rr[2])

            u_ll = prim2cons(SVector(rho_ll_rr[1], v_ll_rr[1], p_ll_rr[1], a_ll_rr[1]),
                             equations)
            u_rr = prim2cons(SVector(rho_ll_rr[2], v_ll_rr[2], p_ll_rr[2], a_ll_rr[2]),
                             equations)

            @test max_abs_speed_naive(u_ll, u_rr, 1, equations) ≈
                  max_abs_speed(u_ll, u_rr, 1, equations)

            @test u_ll ≈ entropy2cons(cons2entropy(u_ll, equations), equations)
        end
    end

    @timed_testset "IdealGlmMhdEquations1D" begin
        for gamma in [1.4, 5 / 3, 7 / 5]
            equations = IdealGlmMhdEquations1D(gamma)

            rho = 42.0
            v1_ll_rr = SVector(0.1, 0.2)
            v2_ll_rr = SVector(0.2, 0.1)
            v3 = 0.0
            p = 0.4
            B1 = 1.01
            B2 = -0.3
            B3 = 0.5

            u_ll = prim2cons(SVector(rho, v1_ll_rr[1], v2_ll_rr[1], v3, p, B1, B2, B3),
                             equations)

            u_rr = prim2cons(SVector(rho, v1_ll_rr[2], v2_ll_rr[2], v3, p, B1, B2, B3),
                             equations)

            @test max_abs_speed_naive(u_ll, u_rr, 1, equations) ≈
                  max_abs_speed(u_ll, u_rr, 1, equations)
        end
    end

    @timed_testset "IdealGlmMhdEquations2D" begin
        for gamma in [1.4, 5 / 3, 7 / 5]
            equations = IdealGlmMhdEquations2D(gamma)

            rho = 42.0
            v1_ll_rr = SVector(0.1, 0.2)
            v2_ll_rr = SVector(0.2, 0.1)
            v3 = 0.0
            p = 0.4
            B1 = 1.01
            B2 = -0.3
            B3 = 0.5
            psi = 0.1

            u_ll = prim2cons(SVector(rho, v1_ll_rr[1], v2_ll_rr[1], v3, p, B1, B2, B3,
                                     psi), equations)

            u_rr = prim2cons(SVector(rho, v1_ll_rr[2], v2_ll_rr[2], v3, p, B1, B2, B3,
                                     psi), equations)

            for orientation in [1, 2]
                @test max_abs_speed_naive(u_ll, u_rr, orientation, equations) ≈
                      max_abs_speed(u_ll, u_rr, 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 max_abs_speed_naive(u_ll, u_rr, normal_direction, equations) ≈
                      max_abs_speed(u_ll, u_rr, normal_direction, equations)
            end
        end
    end

    @timed_testset "IdealGlmMhdEquations3D" begin
        for gamma in [1.4, 5 / 3, 7 / 5]
            equations = IdealGlmMhdEquations3D(gamma)

            rho = 42.0
            v1 = 0.0
            v2_ll_rr = SVector(0.2, 0.1)
            v3_ll_rr = SVector(0.1, 0.2)
            p = 0.4
            B1 = 1.01
            B2 = -0.3
            B3 = 0.5
            psi = 0.1

            u_ll = prim2cons(SVector(rho, v1, v2_ll_rr[1], v3_ll_rr[1], p, B1, B2, B3,
                                     psi), equations)

            u_rr = prim2cons(SVector(rho, v1, v2_ll_rr[2], v3_ll_rr[2], p, B1, B2, B3,
                                     psi), equations)

            for orientation in [1, 2, 3]
                @test max_abs_speed_naive(u_ll, u_rr, orientation, equations) ≈
                      max_abs_speed(u_ll, u_rr, orientation, equations)
            end

            normal_directions = [SVector(1.0, 0.0, 0.0),
                SVector(0.0, 1.0, 0.0),
                SVector(0.0, 0.0, 1.0),
                SVector(0.5, -0.5, 0.2),
                SVector(-1.2, 0.3, 1.4)]

            for normal_direction in normal_directions
                @test max_abs_speed_naive(u_ll, u_rr, normal_direction, equations) ≈
                      max_abs_speed(u_ll, u_rr, normal_direction, equations)
            end
        end
    end

    @timed_testset "IdealGlmMhdMulticomponentEquations1D" begin
        for gamma in [1.4, 5 / 3, 7 / 5]
            equations = IdealGlmMhdMulticomponentEquations1D(gammas = (gamma,
                                                                       gamma),
                                                             gas_constants = (0.5,
                                                                              0.4))

            rho1_ll_rr = SVector(2.0, 1.0)
            rho2_ll_rr = SVector(2.0, 1.0)
            v1_ll_rr = SVector(0.1, 0.2)
            v2_ll_rr = SVector(0.2, 0.1)
            v3 = 0.0
            p = 0.4
            B1 = 1.01
            B2 = -0.3
            B3 = 0.5

            u_ll = prim2cons(SVector(v1_ll_rr[1], v2_ll_rr[1], v3, p, B1, B2, B3,
                                     rho1_ll_rr[1], rho2_ll_rr[1]), equations)
            u_rr = prim2cons(SVector(v1_ll_rr[2], v2_ll_rr[2], v3, p, B1, B2, B3,
                                     rho1_ll_rr[2], rho2_ll_rr[2]), equations)

            @test max_abs_speed_naive(u_ll, u_rr, 1, equations) ≈
                  max_abs_speed(u_ll, u_rr, 1, equations)
        end
    end

    @timed_testset "IdealGlmMhdMulticomponentEquations2D" begin
        for gamma in [1.4, 5 / 3, 7 / 5]
            equations = IdealGlmMhdMulticomponentEquations2D(gammas = (gamma,
                                                                       gamma),
                                                             gas_constants = (0.5,
                                                                              0.4))

            rho1_ll_rr = SVector(0.5, 0.5)
            rho2_ll_rr = SVector(0.5, 0.5)
            v1_ll_rr = SVector(0.1, 0.2)
            v2_ll_rr = SVector(0.2, 0.1)
            v3 = 0.0
            p = 0.4
            B1 = 1.1
            B2 = -0.3
            B3 = 0.4
            psi = 0.1

            u_ll = prim2cons(SVector(v1_ll_rr[1], v2_ll_rr[1], v3, p, B1, B2, B3, psi,
                                     rho1_ll_rr[1], rho2_ll_rr[1]), equations)
            u_rr = prim2cons(SVector(v1_ll_rr[2], v2_ll_rr[2], v3, p, B1, B2, B3, psi,
                                     rho1_ll_rr[2], rho2_ll_rr[2]), equations)

            for orientation in [1, 2]
                @test max_abs_speed_naive(u_ll, u_rr, orientation, equations) ≈
                      max_abs_speed(u_ll, u_rr, orientation, equations)
            end
        end
    end

    @timed_testset "IdealGlmMhdMultiIonEquations2D" begin
        equations = IdealGlmMhdMultiIonEquations2D(gammas = (1.4, 1.667),
                                                   charge_to_mass = (1.0, 2.0))

        B1 = 1.1
        B2 = -0.3
        B3 = 0.4
        rho1_ll_rr = SVector(0.5, 0.5)
        rho2_ll_rr = SVector(0.5, 0.5)
        vx1_ll_rr = SVector(0.1, 0.1)
        vx2_ll_rr = SVector(0.2, 0.2)
        vy1_ll_rr = SVector(0.3, 0.3)
        vy2_ll_rr = SVector(0.4, 0.4)
        vz1 = 0.0
        vz2 = 0.0
        p1 = 0.4
        p2 = 0.4
        psi = 0.1

        u_ll = prim2cons(SVector(B1, B2, B3, rho1_ll_rr[1], rho2_ll_rr[1],
                                 vx1_ll_rr[1], vy1_ll_rr[1], vx2_ll_rr[1], vy2_ll_rr[1],
                                 vz1, vz2, p1, p2, psi), equations)

        u_rr = prim2cons(SVector(B1, B2, B3, rho1_ll_rr[2], rho2_ll_rr[2],
                                 vx1_ll_rr[2], vy1_ll_rr[2], vx2_ll_rr[2], vy2_ll_rr[2],
                                 vz1, vz2, p1, p2, psi), equations)

        for orientation in [1, 2]
            @test max_abs_speed_naive(u_ll, u_rr, orientation, equations) ≈
                  max_abs_speed(u_ll, u_rr, orientation, equations)
        end

        normal_directions = [SVector(1.0, 0.0), SVector(0.0, 1.0), SVector(0.5, -0.5)]
        for normal_direction in normal_directions
            @test max_abs_speed_naive(u_ll, u_rr, normal_direction, equations) ≈
                  max_abs_speed(u_ll, u_rr, normal_direction, equations)
        end
    end

    @timed_testset "IdealGlmMhdMultiIonEquations3D" begin
        equations = IdealGlmMhdMultiIonEquations3D(gammas = (1.4, 1.667),
                                                   charge_to_mass = (1.0, 2.0))

        B1 = 1.1
        B2 = -0.3
        B3 = 0.4
        rho1_ll_rr = SVector(0.5, 0.5)
        rho2_ll_rr = SVector(0.5, 0.5)
        vx1_ll_rr = SVector(0.1, 0.1)
        vx2_ll_rr = SVector(0.2, 0.2)
        vy1_ll_rr = SVector(0.3, 0.3)
        vy2_ll_rr = SVector(0.4, 0.4)
        vz1_ll_rr = SVector(0.5, 0.5)
        vz2_ll_rr = SVector(0.6, 0.6)
        p1 = 0.4
        p2 = 0.4
        psi = 0.1

        u_ll = prim2cons(SVector(B1, B2, B3, rho1_ll_rr[1], rho2_ll_rr[1],
                                 vx1_ll_rr[1], vy1_ll_rr[1], vx2_ll_rr[1], vy2_ll_rr[1],
                                 vz1_ll_rr[1], vz2_ll_rr[1], p1, p2, psi), equations)

        u_rr = prim2cons(SVector(B1, B2, B3, rho1_ll_rr[2], rho2_ll_rr[2],
                                 vx1_ll_rr[2], vy1_ll_rr[2], vx2_ll_rr[2], vy2_ll_rr[2],
                                 vz1_ll_rr[2], vz2_ll_rr[2], p1, p2, psi), equations)

        for orientation in [1, 2, 3]
            @test max_abs_speed_naive(u_ll, u_rr, orientation, equations) ≈
                  max_abs_speed(u_ll, u_rr, orientation, equations)
        end

        normal_directions = [SVector(1.0, 0.0, 0.0),
            SVector(0.0, 1.0, 0.0),
            SVector(0.0, 0.0, 1.0),
            SVector(0.5, -0.5, 0.2),
            SVector(-1.2, 0.3, 1.4)]

        for normal_direction in normal_directions
            @test max_abs_speed_naive(u_ll, u_rr, normal_direction, equations) ≈
                  max_abs_speed(u_ll, u_rr, normal_direction, equations)
        end
    end

    @timed_testset "PolytropicEulerEquations2D" begin
        for gamma in [1.4, 5 / 3, 7 / 5]
            equations = PolytropicEulerEquations2D(gamma, gamma * 0.72)

            rho_ll_rr = SVector(2.0, 2.0)
            v1_ll_rr = SVector(0.1, 0.2)
            v2_ll_rr = SVector(0.4, 0.3)

            u_ll = prim2cons(SVector(rho_ll_rr[1], v1_ll_rr[1], v2_ll_rr[1]), equations)
            u_rr = prim2cons(SVector(rho_ll_rr[2], v1_ll_rr[2], v2_ll_rr[2]), equations)

            for orientation in [1, 2]
                @test max_abs_speed_naive(u_ll, u_rr, orientation, equations) ≈
                      max_abs_speed(u_ll, u_rr, 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 max_abs_speed_naive(u_ll, u_rr, normal_direction, equations) ≈
                      max_abs_speed(u_ll, u_rr, normal_direction, equations)
            end
        end
    end
end

@testitem "Unit: SimpleKronecker" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    N = 3

    NDIMS = 2
    r, s = StartUpDG.nodes(Quad(), N)
    V = StartUpDG.vandermonde(Quad(), N, r, s)
    r1D = StartUpDG.nodes(Line(), N)
    V1D = StartUpDG.vandermonde(Line(), N, r1D)

    x = r + s
    V_kron = Trixi.SimpleKronecker(NDIMS, V1D, eltype(x))

    b = similar(x)
    b_kron = similar(x)
    Trixi.mul!(b, V, x)
    Trixi.mul!(b_kron, V_kron, x)
    @test b ≈ b_kron
end

@testitem "Unit: SummationByPartsOperators + StartUpDG" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    global D = derivative_operator(SummationByPartsOperators.MattssonNordström2004(),
                                   derivative_order = 1,
                                   accuracy_order = 4,
                                   xmin = 0.0, xmax = 1.0,
                                   N = 10)
    dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = D)

    @test StartUpDG.inverse_trace_constant(dg.basis) ≈ 50.8235294117647
end

@testitem "Unit: 1D non-periodic DGMultiMesh" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    # checks whether or not boundary faces are initialized correctly for DGMultiMesh in 1D
    dg = DGMulti(polydeg = 1, element_type = Line(), approximation_type = Polynomial(),
                 surface_integral = SurfaceIntegralWeakForm(flux_central),
                 volume_integral = VolumeIntegralFluxDifferencing(flux_central))
    cells_per_dimension = (1,)
    mesh = DGMultiMesh(dg, cells_per_dimension, periodicity = false)

    @test mesh.boundary_faces[:entire_boundary] == [1, 2]
end

@testitem "Unit: PERK Single p2 Constructors" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    path_coeff_file = mktempdir()
    Trixi.download("https://gist.githubusercontent.com/DanielDoehring/8db0808b6f80e59420c8632c0d8e2901/raw/39aacf3c737cd642636dd78592dbdfe4cb9499af/MonCoeffsS6p2.txt",
                   joinpath(path_coeff_file, "gamma_6.txt"))

    ode_algorithm = Trixi.PairedExplicitRK2(6, path_coeff_file)

    @test isapprox(transpose(ode_algorithm.a_matrix),
                   [0.12405417889682908 0.07594582110317093
                    0.16178873711001726 0.13821126288998273
                    0.16692313960864164 0.2330768603913584
                    0.12281292901258256 0.37718707098741744], atol = 1e-13)

    Trixi.download("https://gist.githubusercontent.com/DanielDoehring/c7a89eaaa857e87dde055f78eae9b94a/raw/2937f8872ffdc08e0dcf444ee35f9ebfe18735b0/Spectrum_2D_IsentropicVortex_CEE.txt",
                   joinpath(path_coeff_file, "spectrum_2d.txt"))

    eig_vals = readdlm(joinpath(path_coeff_file, "spectrum_2d.txt"), ComplexF64)
    tspan = (0.0, 1.0)
    ode_algorithm = Trixi.PairedExplicitRK2(12, tspan, vec(eig_vals))

    @test isapprox(transpose(ode_algorithm.a_matrix),
                   [0.06453812656711647 0.02637096434197444
                    0.09470601372274887 0.041657622640887494
                    0.12332877820069793 0.058489403617483886
                    0.14987015032771522 0.07740257694501203
                    0.1734211495362651 0.0993061231910076
                    0.19261978147948638 0.1255620367023318
                    0.20523340226247055 0.1584029613738931
                    0.20734890429023528 0.20174200480067384
                    0.1913514234997008 0.26319403104575373
                    0.13942836392866081 0.3605716360713392], atol = 1e-13)
end

@testitem "Unit: PERK Single p3 Constructors" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    path_coeff_file = mktempdir()
    Trixi.download("https://gist.githubusercontent.com/warisa-r/0796db36abcd5abe735ac7eebf41b973/raw/32889062fd5dcf7f450748f4f5f0797c8155a18d/a_8_8.txt",
                   joinpath(path_coeff_file, "a_8.txt"))

    ode_algorithm = Trixi.PairedExplicitRK3(8, path_coeff_file)

    @test isapprox(transpose(ode_algorithm.a_matrix),
                   [0.33551678438002486 0.06448322158043965
                    0.49653494442225443 0.10346507941960345
                    0.6496890912144586 0.15031092070647037
                    0.789172498521197 0.21082750147880308
                    0.7522972036571336 0.2477027963428664
                    0.31192569908571666 0.18807430091428337], atol = 1e-13)

    Trixi.download("https://gist.githubusercontent.com/warisa-r/8d93f6a3ae0635e13b9f51ee32ab7fff/raw/54dc5b14be9288e186b745facb5bbcb04d1476f8/EigenvalueList_Refined2.txt",
                   joinpath(path_coeff_file, "spectrum.txt"))

    eig_vals = readdlm(joinpath(path_coeff_file, "spectrum.txt"), ComplexF64)
    tspan = (0.0, 1.0)
    ode_algorithm = Trixi.PairedExplicitRK3(13, tspan, vec(eig_vals))

    @test isapprox(transpose(ode_algorithm.a_matrix),
                   [0.19121164778938382 0.008788355190848427
                    0.28723462747227385 0.012765384448655121
                    0.38017717196008227 0.019822834000382223
                    0.4706748928843403 0.029325107115659724
                    0.557574833668358 0.04242519017349991
                    0.6390917512034328 0.06090823687563831
                    0.7124876770174374 0.08751233490349149
                    0.7736369992226316 0.12636297693551043
                    0.8161315324169078 0.1838684675830921
                    0.7532704453316061 0.2467295546683939
                    0.31168238866709846 0.18831761133290154], atol = 1e-8)
end

@testitem "Unit: PERK Single p4 Constructors" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    path_coeff_file = mktempdir()
    Trixi.download("https://gist.githubusercontent.com/warisa-r/8d93f6a3ae0635e13b9f51ee32ab7fff/raw/54dc5b14be9288e186b745facb5bbcb04d1476f8/EigenvalueList_Refined2.txt",
                   joinpath(path_coeff_file, "spectrum.txt"))

    eig_vals = readdlm(joinpath(path_coeff_file, "spectrum.txt"), ComplexF64)
    tspan = (0.0, 1.0)
    ode_algorithm = Trixi.PairedExplicitRK4(14, tspan, vec(eig_vals))

    @test isapprox(transpose(ode_algorithm.a_matrix),
                   [0.9935765040401348 0.0064234959598652
                    0.9849926812139576 0.0150073187860425
                    0.9731978940975923 0.0268021059024077
                    0.9564664284695985 0.0435335715304015
                    0.9319632992510594 0.0680367007489407
                    0.8955171743167522 0.1044828256832478
                    0.8443975130657495 0.1556024869342504
                    0.7922561745278265 0.2077438254721735
                    0.7722324105428290 0.2277675894571710], atol = 1e-13)
end

@testitem "Unit: Sutherlands Law" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    function mu(u, equations)
        T_ref = 291.15

        R_specific_air = 287.052874
        T = R_specific_air * Trixi.temperature(u, equations)

        C_air = 120.0
        mu_ref_air = 1.827e-5

        return mu_ref_air * (T_ref + C_air) / (T + C_air) * (T / T_ref)^1.5
    end

    function mu_control(u, equations, T_ref, R_specific, C, mu_ref)
        T = R_specific * Trixi.temperature(u, equations)

        return mu_ref * (T_ref + C) / (T + C) * (T / T_ref)^1.5
    end

    # Dry air (values from Wikipedia: https://de.wikipedia.org/wiki/Sutherland-Modell)
    T_ref = 291.15
    C = 120.0 # Sutherland's constant
    R_specific = 287.052874
    mu_ref = 1.827e-5
    prandtl_number() = 0.72
    gamma = 1.4

    equations = CompressibleEulerEquations2D(gamma)
    equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu,
                                                              Prandtl = prandtl_number())
    @test equations_parabolic.gamma == gamma
    @test :gamma in @inferred(propertynames(equations_parabolic))

    # Flow at rest
    u = prim2cons(SVector(1.0, 0.0, 0.0, 1.0), equations_parabolic)

    # Comparison value from https://www.engineeringtoolbox.com/air-absolute-kinematic-viscosity-d_601.html at 18°C
    @test isapprox(mu_control(u, equations_parabolic, T_ref, R_specific, C, mu_ref),
                   1.803e-5, atol = 5e-8)
end

@testitem "Unit: Slope Limiters" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    sl = 1.0
    sr = -1.0

    # Test for code coverage
    dummy = 42
    @test reconstruction_constant(dummy, sl, sr, dummy, dummy, dummy, dummy, dummy) ==
          (sl, sr)

    @test minmod(sl, sr) == 0.0
    @test monotonized_central(sl, sr) == 0.0
    @test superbee(sl, sr) == 0.0
    @test vanleer(sl, sr) == 0.0

    sr = 0.5
    @test minmod(sl, sr) == 0.5
    @test monotonized_central(sl, sr) == 0.75
    @test superbee(sl, sr) == 1.0
    @test isapprox(vanleer(sl, sr), 2 / 3)

    sl = -1.0
    sr = 0.0
    @test minmod(sl, sr) == 0.0
    @test monotonized_central(sl, sr) == 0.0
    @test superbee(sl, sr) == 0.0
    @test vanleer(sl, sr) == 0.0

    sr = -0.8
    @test minmod(sl, sr) == -0.8
    @test monotonized_central(sl, sr) == -0.9
    @test superbee(sl, sr) == -1.0
    @test isapprox(vanleer(sl, sr), -8 / 9)

    # Test symmetry
    @test minmod(sr, sl) == -0.8
    @test monotonized_central(sr, sl) == -0.9
    @test superbee(sr, sl) == -1.0
    @test isapprox(vanleer(sr, sl), -8 / 9)

    sl = 1.0
    sr = 0.0
    @test minmod(sl, sr) == 0.0
    @test monotonized_central(sl, sr) == 0.0
    @test superbee(sl, sr) == 0.0
    @test vanleer(sl, sr) == 0.0

    @test central_slope(sl, sr) == 0.5

    # Test van Leer zero case
    @test vanleer(0.0, 0.0) == 0.0

    sl = -1.0
    sr = -2.0
    @test koren(sl, sr) == -5 / 3
    @test koren(sl, sr) == koren_flipped(sr, sl)
    @test koren_symmetric(sl, sr) == -4 / 3

    sl = 0.0
    @test koren(sl, sr) == 0.0
    @test koren(sl, sr) == koren_flipped(sr, sl)
    @test koren_symmetric(sl, sr) == 0.0

    sr = 2.0
    @test koren(sl, sr) == 0.0
    @test koren(sl, sr) == koren_flipped(sr, sl)
    @test koren_symmetric(sl, sr) == 0.0

    sl = 1.0
    @test koren(sl, sr) == 5 / 3
    @test koren(sl, sr) == koren_flipped(sr, sl)
    @test koren_symmetric(sl, sr) == 4 / 3
end

# Velocity functions are present in many equations and are tested here
@testitem "Unit: Velocity functions for different equations" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    gamma = 1.4
    rho = pi * pi
    pres = sqrt(pi)
    v1, v2, v3 = pi, exp(1.0), exp(pi) # use pi, exp to test with non-trivial numbers
    v_vector = SVector(v1, v2, v3)
    normal_direction_2d = SVector(pi^2, pi^3)
    normal_direction_3d = SVector(normal_direction_2d..., pi^4)
    v_normal_1d = v1 * normal_direction_2d[1]
    v_normal_2d = v1 * normal_direction_2d[1] + v2 * normal_direction_2d[2]
    v_normal_3d = v_normal_2d + v3 * normal_direction_3d[3]

    equations_euler_1d = CompressibleEulerEquations1D(gamma)
    u = prim2cons(SVector(rho, v1, pres), equations_euler_1d)
    @test isapprox(velocity(u, equations_euler_1d), v1)
    orientation = 1 # 1D only has one orientation
    @test isapprox(velocity(u, orientation, equations_euler_1d), v1)

    equations_euler_2d = CompressibleEulerEquations2D(gamma)
    u = prim2cons(SVector(rho, v1, v2, pres), equations_euler_2d)
    @test isapprox(velocity(u, equations_euler_2d), SVector(v1, v2))
    @test isapprox(velocity(u, normal_direction_2d, equations_euler_2d), v_normal_2d)
    for orientation in 1:2
        @test isapprox(velocity(u, orientation, equations_euler_2d),
                       v_vector[orientation])
    end

    equations_euler_3d = CompressibleEulerEquations3D(gamma)
    u = prim2cons(SVector(rho, v1, v2, v3, pres), equations_euler_3d)
    @test isapprox(velocity(u, equations_euler_3d), SVector(v1, v2, v3))
    @test isapprox(velocity(u, normal_direction_3d, equations_euler_3d), v_normal_3d)
    for orientation in 1:3
        @test isapprox(velocity(u, orientation, equations_euler_3d),
                       v_vector[orientation])
    end

    rho1, rho2 = rho, rho * pi # use pi to test with non-trivial numbers
    gammas = (gamma, exp(gamma))
    gas_constants = (0.387, 1.678) # Standard numbers + 0.1

    equations_multi_euler_1d = CompressibleEulerMulticomponentEquations1D(; gammas,
                                                                          gas_constants)
    u = prim2cons(SVector(v1, pres, rho1, rho2), equations_multi_euler_1d)
    @test isapprox(velocity(u, equations_multi_euler_1d), v1)

    equations_multi_euler_2d = CompressibleEulerMulticomponentEquations2D(; gammas,
                                                                          gas_constants)
    u = prim2cons(SVector(v1, v2, pres, rho1, rho2), equations_multi_euler_2d)
    @test isapprox(velocity(u, equations_multi_euler_2d), SVector(v1, v2))
    @test isapprox(velocity(u, normal_direction_2d, equations_multi_euler_2d),
                   v_normal_2d)
    for orientation in 1:2
        @test isapprox(velocity(u, orientation, equations_multi_euler_2d),
                       v_vector[orientation])
    end

    kappa = 0.1 * pi # pi for non-trivial test
    equations_polytropic = PolytropicEulerEquations2D(gamma, kappa)
    u = prim2cons(SVector(rho, v1, v2), equations_polytropic)
    @test isapprox(velocity(u, equations_polytropic), SVector(v1, v2))
    equations_polytropic = CompressibleEulerMulticomponentEquations2D(; gammas,
                                                                      gas_constants)
    u = prim2cons(SVector(v1, v2, pres, rho1, rho2), equations_polytropic)
    @test isapprox(velocity(u, equations_polytropic), SVector(v1, v2))
    @test isapprox(velocity(u, normal_direction_2d, equations_polytropic), v_normal_2d)
    for orientation in 1:2
        @test isapprox(velocity(u, orientation, equations_polytropic),
                       v_vector[orientation])
    end

    B1, B2, B3 = pi^3, pi^4, pi^5
    equations_ideal_mhd_1d = IdealGlmMhdEquations1D(gamma)
    u = prim2cons(SVector(rho, v1, v2, v3, pres, B1, B2, B3), equations_ideal_mhd_1d)
    @test isapprox(velocity(u, equations_ideal_mhd_1d), SVector(v1, v2, v3))
    for orientation in 1:3
        @test isapprox(velocity(u, orientation, equations_ideal_mhd_1d),
                       v_vector[orientation])
    end

    psi = exp(0.1)
    equations_ideal_mhd_2d = IdealGlmMhdEquations2D(gamma)
    u = prim2cons(SVector(rho, v1, v2, v3, pres, B1, B2, B3, psi),
                  equations_ideal_mhd_2d)
    @test isapprox(velocity(u, equations_ideal_mhd_2d), SVector(v1, v2, v3))
    @test isapprox(velocity(u, normal_direction_2d, equations_ideal_mhd_2d),
                   v_normal_2d)
    for orientation in 1:3
        @test isapprox(velocity(u, orientation, equations_ideal_mhd_2d),
                       v_vector[orientation])
    end

    equations_ideal_mhd_3d = IdealGlmMhdEquations3D(gamma)
    u = prim2cons(SVector(rho, v1, v2, v3, pres, B1, B2, B3, psi),
                  equations_ideal_mhd_3d)
    @test isapprox(velocity(u, equations_ideal_mhd_3d), SVector(v1, v2, v3))
    @test isapprox(velocity(u, normal_direction_3d, equations_ideal_mhd_3d),
                   v_normal_3d)
    for orientation in 1:3
        @test isapprox(velocity(u, orientation, equations_ideal_mhd_3d),
                       v_vector[orientation])
    end
end

@testitem "Unit: Pretty_form output for lake_at_rest_error" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    @test Trixi.pretty_form_utf(lake_at_rest_error) == "∑|H₀-(h+b)|"
    @test Trixi.pretty_form_ascii(lake_at_rest_error) == "|H0-(h+b)|"
end

# Ensure consistency for nonconservative fluxes used in the subcell-limiting. Specifically, test
# that flux_noncons_local_structured = flux_noncons_local * flux_noncons_structured.
@testitem "Unit: Nonconservative fluxes for subcell-limiting" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    equations = IdealGlmMhdEquations2D(1.4)
    u_ll = SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1, 0.0)
    u_rr = SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2)

    ## Tests for flux_nonconservative_powell_local_symmetric
    # Implementation for meshes with orientation
    for orientation in 1:2
        flux_noncons = zero(u_ll)
        for noncons in 1:Trixi.n_nonconservative_terms(flux_nonconservative_powell_local_symmetric)
            flux_noncons += flux_nonconservative_powell_local_symmetric(u_ll, 1,
                                                                        equations,
                                                                        Trixi.NonConservativeLocal(),
                                                                        noncons) .*
                            flux_nonconservative_powell_local_symmetric(u_ll, u_rr, 1,
                                                                        equations,
                                                                        Trixi.NonConservativeSymmetric(),
                                                                        noncons)
        end

        @test flux_noncons ≈
              flux_nonconservative_powell_local_symmetric(u_ll, u_rr, 1, equations)
    end

    # Implementation for meshes with normal_direction
    for (orientation, normal_direction) in enumerate((SVector(1.0, 0.0),
                                                      SVector(0.0, 1.0)))
        flux_noncons = zero(u_ll)
        for noncons in 1:Trixi.n_nonconservative_terms(flux_nonconservative_powell_local_symmetric)
            flux_noncons += flux_nonconservative_powell_local_symmetric(u_ll,
                                                                        normal_direction,
                                                                        equations,
                                                                        Trixi.NonConservativeLocal(),
                                                                        noncons) .*
                            flux_nonconservative_powell_local_symmetric(u_ll, u_rr,
                                                                        normal_direction,
                                                                        equations,
                                                                        Trixi.NonConservativeSymmetric(),
                                                                        noncons)
        end

        @test flux_noncons ≈
              flux_nonconservative_powell_local_symmetric(u_ll, u_rr, normal_direction,
                                                          equations)
        @test flux_noncons ≈
              flux_nonconservative_powell_local_symmetric(u_ll, u_rr, orientation,
                                                          equations)
    end

    ## Tests for flux_nonconservative_powell_local_jump
    # Implementation for meshes with orientation
    for orientation in 1:2
        flux_noncons = zero(u_ll)
        for noncons in 1:Trixi.n_nonconservative_terms(flux_nonconservative_powell_local_jump)
            flux_noncons += flux_nonconservative_powell_local_jump(u_ll, 1, equations,
                                                                   Trixi.NonConservativeLocal(),
                                                                   noncons) .*
                            flux_nonconservative_powell_local_jump(u_ll, u_rr, 1,
                                                                   equations,
                                                                   Trixi.NonConservativeJump(),
                                                                   noncons)
        end

        @test flux_noncons ≈
              flux_nonconservative_powell_local_jump(u_ll, u_rr, 1, equations)
    end

    # Implementation for meshes with normal_direction
    for (orientation, normal_direction) in enumerate((SVector(1.0, 0.0),
                                                      SVector(0.0, 1.0)))
        flux_noncons = zero(u_ll)
        for noncons in 1:Trixi.n_nonconservative_terms(flux_nonconservative_powell_local_jump)
            flux_noncons += flux_nonconservative_powell_local_jump(u_ll,
                                                                   normal_direction,
                                                                   equations,
                                                                   Trixi.NonConservativeLocal(),
                                                                   noncons) .*
                            flux_nonconservative_powell_local_jump(u_ll, u_rr,
                                                                   normal_direction,
                                                                   normal_direction,
                                                                   equations,
                                                                   Trixi.NonConservativeJump(),
                                                                   noncons)
        end

        @test flux_noncons ≈
              flux_nonconservative_powell_local_jump(u_ll, u_rr, normal_direction,
                                                     equations)
        @test flux_noncons ≈
              flux_nonconservative_powell_local_jump(u_ll, u_rr, orientation,
                                                     equations)
    end
end

@testitem "Unit: SparseConnectivityTracer FiniteDiff Jacobian" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    ###############################################################################
    ### equations, solver, mesh ###

    advection_velocities = (0.2, -0.7)
    equations = LinearScalarAdvectionEquation2D(advection_velocities)

    float_type = Float64 # Datatype for the actual simulation
    solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, RealT = float_type)

    coordinates_min = (-1.0, -1.0)
    coordinates_max = (1.0, 1.0)

    mesh = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level = 4,
                    periodicity = true)
    ###############################################################################
    ### semidiscretization for sparsity detection ###

    jac_detector = TracerSparsityDetector()
    # We need to construct the semidiscretization with the correct
    # sparsity-detection ready datatype, which is retrieved here
    jac_eltype = jacobian_eltype(float_type, jac_detector)

    # Semidiscretization for sparsity pattern detection
    semi_jac_type = SemidiscretizationHyperbolic(mesh, equations,
                                                 initial_condition_convergence_test,
                                                 solver;
                                                 boundary_conditions = boundary_condition_periodic,
                                                 uEltype = jac_eltype) # Need to supply Jacobian element type

    tspan = (0.0, 1.0) # Re-used for wrapping `rhs` below

    # Call `semidiscretize` to create the ODE problem to have access to the
    # initial condition based on which the sparsity pattern is computed
    ode_jac_type = semidiscretize(semi_jac_type, tspan)
    u0_ode = ode_jac_type.u0
    du_ode = similar(u0_ode)

    ###############################################################################
    ### Compute the Jacobian sparsity pattern ###

    # Wrap the `Trixi.rhs_hyperbolic!` function to match the signature `f!(du, u)`, see
    # https://adrianhill.de/SparseConnectivityTracer.jl/stable/user/api/#ADTypes.jacobian_sparsity
    rhs_jac_type! = function (du_ode, u0_ode)
        Trixi.rhs_hyperbolic!(du_ode, u0_ode, semi_jac_type, tspan[1])
    end

    jac_prototype = jacobian_sparsity(rhs_jac_type!, du_ode, u0_ode, jac_detector)

    coloring_prob = ColoringProblem(; structure = :nonsymmetric, partition = :column)
    coloring_alg = GreedyColoringAlgorithm(; decompression = :direct)
    coloring_result = coloring(jac_prototype, coloring_prob, coloring_alg)
    coloring_vec = column_colors(coloring_result)

    ###############################################################################
    ### float-type semidiscretization ###

    semi_float_type = SemidiscretizationHyperbolic(mesh, equations,
                                                   initial_condition_convergence_test,
                                                   solver;
                                                   boundary_conditions = boundary_condition_periodic)

    ode_float_type = semidiscretize(semi_float_type, tspan)
    u0_ode = ode_float_type.u0
    du_ode = similar(u0_ode)
    N = length(u0_ode)

    @test Trixi.default_rhs(semi_float_type) === Trixi.rhs_hyperbolic!

    rhs_float_type! = function (du_ode, u0_ode)
        Trixi.rhs_hyperbolic!(du_ode, u0_ode, semi_float_type, tspan[1])
    end

    ###############################################################################
    ### sparsity-aware finite diff ###

    jac_sparse_finite_diff = spzeros(N, N)
    finite_difference_jacobian!(jac_sparse_finite_diff, rhs_float_type!,
                                u0_ode, sparsity = jac_prototype,
                                colorvec = coloring_vec)

    jac_finite_diff = jacobian_fd(semi_float_type)

    @test isapprox(jac_finite_diff, jac_sparse_finite_diff; rtol = 5e-8)
    @test isapprox(jac_finite_diff, Matrix(jac_sparse_finite_diff); rtol = 5e-8)
    @test isapprox(sparse(jac_finite_diff), jac_sparse_finite_diff; rtol = 5e-8)
end

@testitem "Unit: Parabolic-Hyperbolic Problem Sparsity Pattern" setup=[Setup, UnitTests] tags=[:misc_part1] begin

    # Poor-mans rebuild of `SplitODEProblem` from SciML
    function rhs_hyperbolic_parabolic!(du_ode, u_ode,
                                       semi::SemidiscretizationHyperbolicParabolic, t)
        du_para = similar(du_ode) # This obviously allocates, but fine for this test
        Trixi.rhs_hyperbolic!(du_ode, u_ode, semi, t)
        Trixi.rhs_parabolic!(du_para, u_ode, semi, t)

        Trixi.@threaded for i in eachindex(du_ode)
            du_ode[i] = du_ode[i] + du_para[i]
        end
        return nothing
    end

    ###############################################################################
    ### equations, solver, mesh ###

    advection_velocity = 1.5
    equations_hyperbolic = LinearScalarAdvectionEquation1D(advection_velocity)
    diffusivity() = 5.0e-2
    equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations_hyperbolic)

    solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)

    coordinates_min = -1.0
    coordinates_max = 1.0

    mesh = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level = 4,
                    periodicity = true)

    ###############################################################################
    ### semidiscretization for sparsity detection ###

    jac_detector = TracerSparsityDetector()
    # We need to construct the semidiscretization with the correct
    # sparsity-detection ready datatype, which is retrieved here
    jac_eltype = jacobian_eltype(real(solver), jac_detector)

    # Semidiscretization for sparsity pattern detection
    semi_jac_type = SemidiscretizationHyperbolicParabolic(mesh,
                                                          (equations_hyperbolic,
                                                           equations_parabolic),
                                                          initial_condition_convergence_test,
                                                          solver;
                                                          boundary_conditions = (boundary_condition_periodic,
                                                                                 boundary_condition_periodic),
                                                          uEltype = jac_eltype) # Need to supply Jacobian element type

    @test_throws ArgumentError Trixi.default_rhs(semi_jac_type)

    tspan = (0.0, 1.5) # Re-used for wrapping `rhs` below

    # Call `semidiscretize` to create the ODE problem to have access to the
    # initial condition based on which the sparsity pattern is computed
    ode_jac_type = semidiscretize(semi_jac_type, tspan)
    u0_ode = ode_jac_type.u0
    du_ode = similar(u0_ode)

    ###############################################################################
    ### Compute the Jacobian sparsity pattern ###

    # Only the parabolic part of the `SplitODEProblem` is treated implicitly so we only need the parabolic Jacobian, see
    # https://docs.sciml.ai/DiffEqDocs/stable/types/split_ode_types/#SciMLBase.SplitFunction
    # Thus, we perform sparsity detection on `rhs_parabolic!` only,
    # which is equivalent to doing sparsity detection on the entire hyperbolic-parabolic problem,
    # at least for the DGSEM & Bassi-Rebay 1 parabolic solver.
    # This test validates this.

    # Wrap the `Trixi.rhs_parabolic!` function to match the signature `f!(du, u)`, see
    # https://adrianhill.de/SparseConnectivityTracer.jl/stable/user/api/#ADTypes.jacobian_sparsity
    rhs_parabolic_wrapped! = (du_ode, u0_ode) -> Trixi.rhs_parabolic!(du_ode,
                                                                      u0_ode,
                                                                      semi_jac_type,
                                                                      tspan[1])

    jac_prototype_parabolic = jacobian_sparsity(rhs_parabolic_wrapped!,
                                                du_ode, u0_ode,
                                                jac_detector)

    ###############################################################################
    ### Compare sparsity pattern detected using `rhs_parabolic!` only to ###
    ### sparsity pattern detected on the combined hyperbolic-parabolic RHS ###

    rhs_hyp_para_wrapped! = (du_ode, u0_ode) -> rhs_hyperbolic_parabolic!(du_ode,
                                                                          u0_ode,
                                                                          semi_jac_type,
                                                                          tspan[1])

    jac_prototype_hyperbolic_parabolic = jacobian_sparsity(rhs_hyp_para_wrapped!,
                                                           du_ode, u0_ode,
                                                           jac_detector)

    # Given that the stencil for the BR1 parabolic solver is for the DGSEM always larger than that of a hyperbolic solver,
    # the sparsity pattern of the parabolic part of a hyperbolic-parabolic problem always includes the hyperbolic one
    @test jac_prototype_parabolic == jac_prototype_hyperbolic_parabolic
end

@testitem "Unit: TreeMesh and StructuredMesh boundary condition argument checks" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    solver = DGSEM(polydeg = 1)
    ic = initial_condition_convergence_test
    bc = boundary_condition_periodic
    bc_dn = boundary_condition_do_nothing
    # 1D
    eq1d = LinearScalarAdvectionEquation1D(1.0)
    tree_mesh1d_periodic = TreeMesh((-1.0,), (1.0,), initial_refinement_level = 1,
                                    periodicity = true)
    structured_mesh1d_periodic = StructuredMesh((4,), (-1.0,), (1.0,),
                                                periodicity = true)
    for mesh1d_periodic in (tree_mesh1d_periodic,
                            structured_mesh1d_periodic)
        # Passing Tuples and arrays is not allowed
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh1d_periodic, eq1d,
                                                                ic, solver;
                                                                boundary_conditions = (bc,
                                                                                       bc))
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh1d_periodic, eq1d,
                                                                ic, solver;
                                                                boundary_conditions = [bc,
                                                                    bc])
        @test_nowarn SemidiscretizationHyperbolic(mesh1d_periodic, eq1d,
                                                  initial_condition_convergence_test,
                                                  solver;
                                                  boundary_conditions = bc)
        # Not passing periodic boundary conditions in NamedTuple for periodic mesh is allowed
        @test_nowarn SemidiscretizationHyperbolic(mesh1d_periodic, eq1d,
                                                  ic, solver;
                                                  boundary_conditions = (;
                                                                         x_neg = bc,))
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh1d_periodic, eq1d,
                                                                ic, solver;
                                                                boundary_conditions = (;
                                                                                       x_neg = bc_dn,))
        # Wrong keys NamedTuple
        @test_throws ErrorException SemidiscretizationHyperbolic(mesh1d_periodic, eq1d,
                                                                 ic, solver;
                                                                 boundary_conditions = (;
                                                                                        x_neg = bc,
                                                                                        y_pos = bc))
    end
    # non-periodic mesh
    tree_mesh1d_nonperiodic = TreeMesh((-1.0,), (1.0,), initial_refinement_level = 1,
                                       periodicity = false)
    structured_mesh1d_nonperiodic = StructuredMesh((4,), (-1.0,), (1.0,),
                                                   periodicity = false)
    for mesh1d_nonperiodic in (tree_mesh1d_nonperiodic,
                               structured_mesh1d_nonperiodic)
        @test_nowarn SemidiscretizationHyperbolic(mesh1d_nonperiodic,
                                                  eq1d, ic, solver;
                                                  boundary_conditions = (;
                                                                         x_neg = bc_dn,
                                                                         x_pos = bc_dn))
        # periodic boundary conditions for non-periodic mesh is not allowed
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh1d_nonperiodic,
                                                                eq1d, ic, solver;
                                                                boundary_conditions = bc)
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh1d_nonperiodic,
                                                                eq1d, ic, solver;
                                                                boundary_conditions = (;
                                                                                       x_neg = bc_dn,
                                                                                       x_pos = bc))
        # not passing non-periodic boundary conditions for non-periodic mesh is not allowed
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh1d_nonperiodic,
                                                                eq1d, ic, solver;
                                                                boundary_conditions = (;
                                                                                       x_neg = bc_dn,))
    end
    # 2D
    eq2d = LinearScalarAdvectionEquation2D((1.0, -1.0))
    tree_mesh2d_periodic = TreeMesh((-1.0, -1.0), (1.0, 1.0),
                                    initial_refinement_level = 1,
                                    periodicity = true)
    structured_mesh2d_periodic = StructuredMesh((4, 4), (-1.0, -1.0), (1.0, 1.0),
                                                periodicity = true)
    for mesh2d_periodic in (tree_mesh2d_periodic,
                            structured_mesh2d_periodic)
        # Passing Tuples and arrays is not allowed
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh2d_periodic, eq2d,
                                                                ic, solver;
                                                                boundary_conditions = (bc,
                                                                                       bc,
                                                                                       bc,
                                                                                       bc))
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh2d_periodic, eq2d,
                                                                ic, solver;
                                                                boundary_conditions = [bc,
                                                                    bc, bc, bc])
        @test_nowarn SemidiscretizationHyperbolic(mesh2d_periodic, eq2d,
                                                  ic, solver;
                                                  boundary_conditions = bc)
        @test_nowarn SemidiscretizationHyperbolic(mesh2d_periodic, eq2d,
                                                  ic, solver;
                                                  boundary_conditions = (; x_neg = bc,
                                                                         x_pos = bc,
                                                                         y_neg = bc,
                                                                         y_pos = bc))
        # Not passing periodic boundary conditions in NamedTuple for periodic mesh is allowed
        @test_nowarn SemidiscretizationHyperbolic(mesh2d_periodic, eq2d,
                                                  ic, solver;
                                                  boundary_conditions = (;
                                                                         x_neg = bc,))
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh2d_periodic, eq2d,
                                                                ic, solver;
                                                                boundary_conditions = (;
                                                                                       x_neg = bc_dn,))
        # Wrong keys NamedTuple
        @test_throws ErrorException SemidiscretizationHyperbolic(mesh2d_periodic, eq2d,
                                                                 ic, solver;
                                                                 boundary_conditions = (;
                                                                                        x_neg = bc,
                                                                                        x_pos = bc,
                                                                                        z_neg = bc,
                                                                                        z_pos = bc))
    end
    # non-periodic mesh
    tree_mesh2d_nonperiodic = TreeMesh((-1.0, -1.0), (1.0, 1.0),
                                       initial_refinement_level = 1,
                                       periodicity = false)
    structured_mesh2d_nonperiodic = StructuredMesh((4, 4), (-1.0, -1.0), (1.0, 1.0),
                                                   periodicity = false)
    for mesh2d_nonperiodic in (tree_mesh2d_nonperiodic,
                               structured_mesh2d_nonperiodic)
        @test_nowarn SemidiscretizationHyperbolic(mesh2d_nonperiodic,
                                                  eq2d, ic, solver;
                                                  boundary_conditions = (;
                                                                         x_neg = bc_dn,
                                                                         x_pos = bc_dn,
                                                                         y_neg = bc_dn,
                                                                         y_pos = bc_dn))
        # periodic boundary conditions for non-periodic mesh is not allowed
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh2d_nonperiodic,
                                                                eq2d, ic, solver;
                                                                boundary_conditions = bc)
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh2d_nonperiodic,
                                                                eq2d, ic, solver;
                                                                boundary_conditions = (;
                                                                                       x_neg = bc_dn,
                                                                                       x_pos = bc,
                                                                                       y_neg = bc,
                                                                                       y_pos = bc))
        # not passing non-periodic boundary conditions for non-periodic mesh is not allowed
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh2d_nonperiodic,
                                                                eq2d, ic, solver;
                                                                boundary_conditions = (;
                                                                                       x_neg = bc_dn,
                                                                                       x_pos = bc_dn,
                                                                                       y_neg = bc_dn))
    end
    # partially periodic
    tree_mesh2d_partial_periodic = TreeMesh((-1.0, -1.0), (1.0, 1.0),
                                            initial_refinement_level = 1,
                                            periodicity = (true, false))
    structured_mesh2d_partial_periodic = StructuredMesh((4, 4), (-1.0, -1.0),
                                                        (1.0, 1.0),
                                                        periodicity = (true, false))
    for mesh2d_partial_periodic in (tree_mesh2d_partial_periodic,
                                    structured_mesh2d_partial_periodic)
        # Specifying all boundary conditions is allowed
        @test_nowarn SemidiscretizationHyperbolic(mesh2d_partial_periodic,
                                                  eq2d, ic, solver;
                                                  boundary_conditions = (;
                                                                         x_neg = bc,
                                                                         x_pos = bc,
                                                                         y_neg = bc_dn,
                                                                         y_pos = bc_dn))
        # Only specifying non-periodic boundary conditions is allowed when using NamedTuple
        @test_nowarn SemidiscretizationHyperbolic(mesh2d_partial_periodic,
                                                  eq2d, ic, solver;
                                                  boundary_conditions = (;
                                                                         y_neg = bc_dn,
                                                                         y_pos = bc_dn))
        # For partially periodic mesh, need to specify separate boundary conditions
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh2d_partial_periodic,
                                                                eq2d, ic, solver;
                                                                boundary_conditions = bc_dn)
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh2d_partial_periodic,
                                                                eq2d, ic, solver;
                                                                boundary_conditions = bc)
        # Non-periodic boundary condition on periodic direction
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh2d_partial_periodic,
                                                                eq2d, ic, solver;
                                                                boundary_conditions = (;
                                                                                       x_neg = bc_dn,
                                                                                       x_pos = bc_dn,
                                                                                       y_neg = bc_dn,
                                                                                       y_pos = bc_dn))
    end
    # 3D
    eq3d = LinearScalarAdvectionEquation3D((1.0, 1.0, -1.0))
    tree_mesh3d_periodic = TreeMesh((-1.0, -1.0, -1.0), (1.0, 1.0, 1.0),
                                    initial_refinement_level = 1,
                                    periodicity = true)
    structured_mesh3d_periodic = StructuredMesh((4, 4, 4), (-1.0, -1.0, -1.0),
                                                (1.0, 1.0, 1.0),
                                                periodicity = true)
    for mesh3d_periodic in (tree_mesh3d_periodic,
                            structured_mesh3d_periodic)
        # Passing Tuples and arrays is not allowed
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh3d_periodic, eq3d,
                                                                ic, solver;
                                                                boundary_conditions = (bc,
                                                                                       bc,
                                                                                       bc,
                                                                                       bc,
                                                                                       bc,
                                                                                       bc))
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh3d_periodic, eq3d,
                                                                ic, solver;
                                                                boundary_conditions = [bc,
                                                                    bc, bc, bc, bc, bc])
        @test_nowarn SemidiscretizationHyperbolic(mesh3d_periodic, eq3d,
                                                  ic, solver;
                                                  boundary_conditions = bc)
        @test_nowarn SemidiscretizationHyperbolic(mesh3d_periodic, eq3d,
                                                  ic, solver;
                                                  boundary_conditions = (; x_neg = bc,
                                                                         x_pos = bc,
                                                                         y_neg = bc,
                                                                         y_pos = bc,
                                                                         z_neg = bc,
                                                                         z_pos = bc))
        # Not passing periodic boundary conditions in NamedTuple for periodic mesh is allowed
        @test_nowarn SemidiscretizationHyperbolic(mesh3d_periodic, eq3d,
                                                  ic, solver;
                                                  boundary_conditions = (;
                                                                         x_neg = bc,))
        # Passing non-periodic boundary conditions for periodic mesh is not allowed
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh3d_periodic, eq3d,
                                                                ic, solver;
                                                                boundary_conditions = (;
                                                                                       x_neg = bc_dn,))
        # Wrong keys NamedTuple
        @test_throws ErrorException SemidiscretizationHyperbolic(mesh3d_periodic, eq3d,
                                                                 ic, solver;
                                                                 boundary_conditions = (;
                                                                                        x_neg = bc,
                                                                                        x_pos = bc,
                                                                                        y_neg = bc,
                                                                                        y_pos = bc,
                                                                                        z_neg = bc,
                                                                                        pos = bc))
    end
    # non-periodic mesh
    tree_mesh3d_nonperiodic = TreeMesh((-1.0, -1.0, -1.0), (1.0, 1.0, 1.0),
                                       initial_refinement_level = 1,
                                       periodicity = false)
    structured_mesh3d_nonperiodic = StructuredMesh((4, 4, 4), (-1.0, -1.0, -1.0),
                                                   (1.0, 1.0, 1.0),
                                                   periodicity = false)
    for mesh3d_nonperiodic in (tree_mesh3d_nonperiodic,
                               structured_mesh3d_nonperiodic)
        # Passing all non-periodic boundary conditions for non-periodic mesh is allowed
        @test_nowarn SemidiscretizationHyperbolic(mesh3d_nonperiodic,
                                                  eq3d, ic, solver;
                                                  boundary_conditions = (;
                                                                         x_neg = bc_dn,
                                                                         x_pos = bc_dn,
                                                                         y_neg = bc_dn,
                                                                         y_pos = bc_dn,
                                                                         z_neg = bc_dn,
                                                                         z_pos = bc_dn))
        # periodic boundary conditions for non-periodic mesh is not allowed
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh3d_nonperiodic,
                                                                eq3d, ic, solver;
                                                                boundary_conditions = bc)
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh3d_nonperiodic,
                                                                eq3d, ic, solver;
                                                                boundary_conditions = (;
                                                                                       x_neg = bc_dn,
                                                                                       x_pos = bc,
                                                                                       y_neg = bc,
                                                                                       y_pos = bc,
                                                                                       z_neg = bc,
                                                                                       z_pos = bc))
        # not passing non-periodic boundary conditions for non-periodic mesh is not allowed
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh3d_nonperiodic,
                                                                eq3d, ic, solver;
                                                                boundary_conditions = (;
                                                                                       x_neg = bc_dn,
                                                                                       x_pos = bc_dn,
                                                                                       y_neg = bc_dn,
                                                                                       y_pos = bc_dn,
                                                                                       z_neg = bc_dn))
    end
    # partially periodic
    tree_mesh3d_partial_periodic = TreeMesh((-1.0, -1.0, -1.0), (1.0, 1.0, 1.0),
                                            initial_refinement_level = 1,
                                            periodicity = (false, true, true))
    structured_mesh3d_partial_periodic = StructuredMesh((4, 4, 4), (-1.0, -1.0, -1.0),
                                                        (1.0, 1.0, 1.0),
                                                        periodicity = (false, true,
                                                                       true))
    for mesh3d_partial_periodic in (tree_mesh3d_partial_periodic,
                                    structured_mesh3d_partial_periodic)
        # Specifying all boundary conditions is allowed
        @test_nowarn SemidiscretizationHyperbolic(mesh3d_partial_periodic,
                                                  eq3d, ic, solver;
                                                  boundary_conditions = (;
                                                                         x_neg = bc_dn,
                                                                         x_pos = bc_dn,
                                                                         y_neg = bc,
                                                                         y_pos = bc,
                                                                         z_neg = bc,
                                                                         z_pos = bc))
        # Only specifying non-periodic boundary conditions is allowed when using NamedTuple
        @test_nowarn SemidiscretizationHyperbolic(mesh3d_partial_periodic,
                                                  eq3d, ic, solver;
                                                  boundary_conditions = (;
                                                                         x_neg = bc_dn,
                                                                         x_pos = bc_dn))
        # For partially periodic mesh, need to specify separate boundary conditions
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh3d_partial_periodic,
                                                                eq3d, ic, solver;
                                                                boundary_conditions = bc_dn)
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh3d_partial_periodic,
                                                                eq3d, ic, solver;
                                                                boundary_conditions = bc)
        # Non-periodic boundary condition on periodic direction
        @test_throws ArgumentError SemidiscretizationHyperbolic(mesh3d_partial_periodic,
                                                                eq3d, ic, solver;
                                                                boundary_conditions = (;
                                                                                       x_neg = bc_dn,
                                                                                       x_pos = bc_dn,
                                                                                       y_neg = bc_dn,
                                                                                       y_pos = bc_dn,
                                                                                       z_neg = bc,
                                                                                       z_pos = bc_dn))
    end
end

@testitem "Unit: ndims function for SemidiscretizaionHyperbolicSplit" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    gamma = 1004 / 717
    equations = CompressibleEulerEquations2D(gamma)

    polydeg = 2
    basis = LobattoLegendreBasis(polydeg)

    volume_integral_explicit = VolumeIntegralFluxDifferencing(flux_ranocha)
    solver_explicit = DGSEM(basis, flux_ranocha, volume_integral_explicit)

    volume_integral_implicit = VolumeIntegralFluxDifferencing(flux_ranocha)
    solver_implicit = DGSEM(basis, flux_ranocha, volume_integral_implicit)

    coordinates_min = (0.0, 0.0)
    coordinates_max = (20_000.0, 10_000.0)
    trees_per_dimension = (16, 8)
    mesh = P4estMesh(trees_per_dimension; polydeg = polydeg,
                     coordinates_min = coordinates_min,
                     coordinates_max = coordinates_max,
                     periodicity = (true, false), initial_refinement_level = 0)

    boundary_conditions = (; y_neg = boundary_condition_slip_wall,
                           y_pos = boundary_condition_slip_wall)

    initial_condition = initial_condition_convergence_test

    semi = SemidiscretizationHyperbolicSplit(mesh,
                                             (equations, equations),
                                             initial_condition,
                                             (solver_implicit, solver_explicit);
                                             boundary_conditions = (boundary_conditions,
                                                                    boundary_conditions),
                                             source_terms = (nothing, nothing),)

    @test Trixi.ndims(semi) == 2
end

@testitem "Unit: Unified mesh constructor signatures (StructuredMesh)" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    # 1D: keyword interface (2^2 = 4 cells per dimension)
    mesh_1d_ref = StructuredMesh((4,), (-1.0,), (1.0,))
    mesh_1d_kw = StructuredMesh(; coordinates_min = (-1.0,), coordinates_max = (1.0,),
                                refinement_level = 2)
    @test mesh_1d_ref.cells_per_dimension == mesh_1d_kw.cells_per_dimension

    # 2D: keyword interface
    mesh_2d_ref = StructuredMesh((4, 4), (-1.0, -1.0), (1.0, 1.0))
    mesh_2d_kw = StructuredMesh(; coordinates_min = (-1.0, -1.0),
                                coordinates_max = (1.0, 1.0),
                                refinement_level = 2)
    @test mesh_2d_ref.cells_per_dimension == mesh_2d_kw.cells_per_dimension

    # 3D: keyword interface
    mesh_3d_ref = StructuredMesh((4, 4, 4), (-1.0, -1.0, -1.0), (1.0, 1.0, 1.0))
    mesh_3d_kw = StructuredMesh(; coordinates_min = (-1.0, -1.0, -1.0),
                                coordinates_max = (1.0, 1.0, 1.0),
                                refinement_level = 2)
    @test mesh_3d_ref.cells_per_dimension == mesh_3d_kw.cells_per_dimension
    @test_throws ArgumentError StructuredMesh(; coordinates_min = (-1.0, -1.0),
                                              coordinates_max = (1.0, 1.0, 1.0),
                                              refinement_level = 2)
end

@testitem "Unit: Unified mesh constructor signatures (DGMultiMesh)" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    dg_1d = DGMulti(polydeg = 2, element_type = Line(),
                    approximation_type = Polynomial(),
                    surface_integral = SurfaceIntegralWeakForm(flux_central),
                    volume_integral = VolumeIntegralFluxDifferencing(flux_central))

    # 1D: keyword interface (2^2 = 4 elements)
    mesh_1d_ref = DGMultiMesh(dg_1d, (4,))
    mesh_1d_kw = DGMultiMesh(dg_1d; coordinates_min = (-1.0,), coordinates_max = (1.0,),
                             refinement_level = 2)
    @test mesh_1d_ref.md.num_elements == mesh_1d_kw.md.num_elements

    dg_2d = DGMulti(polydeg = 2, element_type = Quad(),
                    approximation_type = Polynomial(),
                    surface_integral = SurfaceIntegralWeakForm(flux_central),
                    volume_integral = VolumeIntegralFluxDifferencing(flux_central))

    # 2D: keyword interface
    mesh_2d_ref = DGMultiMesh(dg_2d, (4, 4))
    mesh_2d_kw = DGMultiMesh(dg_2d; coordinates_min = (-1.0, -1.0),
                             coordinates_max = (1.0, 1.0), refinement_level = 2)
    @test mesh_2d_ref.md.num_elements == mesh_2d_kw.md.num_elements
    @test_throws ArgumentError DGMultiMesh(dg_2d; coordinates_min = (-1.0, -1.0),
                                           coordinates_max = (1.0, 1.0, 1.0),
                                           refinement_level = 2)
end

@testitem "Unit: TreeMesh" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    for NDIMS in 1:3
        coords_min = ntuple(_ -> -1.0, NDIMS)
        coords_max = ntuple(_ -> 1.0, NDIMS)
        mesh = TreeMesh(coords_min, coords_max; initial_refinement_level = 2)
        expected_capacity = sum((2^NDIMS)^l for l in 0:2)
        @test @inferred(Trixi.ncells(mesh)) == 2^(NDIMS * 2)
        @test mesh.tree.capacity == expected_capacity
        @test mesh.tree.capacity >= mesh.tree.length
    end
end

@testitem "Unit: TreeMesh auto-growth matches large-capacity tree" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    for NDIMS in 1:2
        RealT = Float64
        TreeType = Trixi.SerialTree{NDIMS, RealT}
        domain_center = SVector{NDIMS, RealT}(ntuple(_ -> 0.0, NDIMS))
        domain_length = convert(RealT, 2.0)

        # Reference: large capacity, no growth needed
        mesh_ref = TreeMesh{NDIMS, TreeType, RealT}(10_000, domain_center,
                                                    domain_length)
        Trixi.initialize!(mesh_ref, 3, (), ())

        # Small: deliberately tiny initial capacity, must grow
        mesh_small = TreeMesh{NDIMS, TreeType, RealT}(2, domain_center,
                                                      domain_length)
        Trixi.initialize!(mesh_small, 3, (), ())

        # Post-construction AMR: refine all leaf cells once on both trees
        Trixi.refine!(mesh_ref.tree)
        Trixi.refine!(mesh_small.tree)

        tr = mesh_ref.tree
        ts = mesh_small.tree

        @test ts.length == tr.length
        @test ts.capacity >= ts.length
        @test ts.parent_ids[1:(ts.length)] == tr.parent_ids[1:(tr.length)]
        @test ts.child_ids[:, 1:(ts.length)] == tr.child_ids[:, 1:(tr.length)]
        @test ts.neighbor_ids[:, 1:(ts.length)] == tr.neighbor_ids[:, 1:(tr.length)]
        @test ts.levels[1:(ts.length)] == tr.levels[1:(tr.length)]
        @test ts.coordinates[:, 1:(ts.length)] ≈ tr.coordinates[:, 1:(tr.length)]
        @test ts.original_cell_ids[1:(ts.length)] == tr.original_cell_ids[1:(tr.length)]
        # Wrapped matrix sizes must match current capacity
        @test size(ts.child_ids) == (2^NDIMS, ts.capacity + 1)
        @test size(ts.neighbor_ids) == (2 * NDIMS, ts.capacity + 1)
        @test size(ts.coordinates) == (NDIMS, ts.capacity + 1)
    end
end

@testitem "Unit: load_mesh derives TreeMesh capacity from n_cells" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    mktempdir() do dir
        mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0);
                        initial_refinement_level = 2)
        mesh_file = Trixi.save_mesh_file(mesh, dir)

        loaded = Trixi.load_mesh_serial(mesh_file; RealT = Float64)
        @test loaded.tree.capacity == mesh.tree.length
        @test loaded.tree.length == mesh.tree.length
        @test loaded.tree.parent_ids[1:(loaded.tree.length)] ==
              mesh.tree.parent_ids[1:(mesh.tree.length)]
        @test loaded.tree.child_ids[:, 1:(loaded.tree.length)] ==
              mesh.tree.child_ids[:, 1:(mesh.tree.length)]
        @test loaded.tree.neighbor_ids[:, 1:(loaded.tree.length)] ==
              mesh.tree.neighbor_ids[:, 1:(mesh.tree.length)]
        @test loaded.tree.levels[1:(loaded.tree.length)] ==
              mesh.tree.levels[1:(mesh.tree.length)]
        @test loaded.tree.coordinates[:, 1:(loaded.tree.length)] ≈
              mesh.tree.coordinates[:, 1:(mesh.tree.length)]
        @test loaded.tree.center_level_0 == mesh.tree.center_level_0
        @test loaded.tree.length_level_0 == mesh.tree.length_level_0
        @test loaded.tree.periodicity == mesh.tree.periodicity
        @test loaded.current_filename == mesh_file
        @test loaded.unsaved_changes == false
    end
end

@testitem "Unit: removed TreeMesh capacity keyword is rejected" setup=[Setup, UnitTests] tags=[:misc_part1] begin
    removed_kw = Symbol("n_cells", "_max")
    kwargs = (; initial_refinement_level = 1, removed_kw => 10)
    @test_throws MethodError TreeMesh((-1.0,), (1.0,); kwargs...)
end

@testitem "Unit: Euler admissible projection for PositivityPreservingLimiterLiuZhang" setup=[
    Setup,
    UnitTests
] tags=[:misc_part1] begin
    @testset "1D projection with different density and internal energy floors" begin
        equations = CompressibleEulerEquations1D(1.4)
        u = SVector(0.5, 1.0, 0.1)
        lower_bounds = (1.0, 0.1)
        variables = (density, energy_internal)

        u_projected = Trixi.project_to_admissible_set(u, lower_bounds, variables, equations)
        arithmetic_tol = Trixi.euler_arithmetic_tol(lower_bounds[1], lower_bounds[2])

        @test u_projected[1] >= lower_bounds[1]
        @test energy_internal(u_projected, equations) >= lower_bounds[2] - arithmetic_tol

        @test u_projected[1] ≈ 1.0
        @test u_projected[2]≈0.7709169970592479 rtol=1e-12
        @test u_projected[3]≈0.39715650817742415 rtol=1e-12
    end

    @testset "Consistency between 1D and 2D projections when v2 = 0" begin
        equations_1d = CompressibleEulerEquations1D(1.4)
        equations_2d = CompressibleEulerEquations2D(1.4)
        u_1d = SVector(0.5, 1.0, 0.1)
        u_2d = SVector(0.5, 1.0, 0.0, 0.1)
        lower_bounds = (1.0, 0.1)
        variables = (density, energy_internal)

        u_projected_1d = Trixi.project_to_admissible_set(u_1d, lower_bounds, variables,
                                                         equations_1d)
        u_projected_2d = Trixi.project_to_admissible_set(u_2d, lower_bounds, variables,
                                                         equations_2d)

        @test u_projected_2d[1] ≈ u_projected_1d[1]
        @test u_projected_2d[2] ≈ u_projected_1d[2]
        @test u_projected_2d[4] ≈ u_projected_1d[3]
        @test u_projected_2d[3] == 0.0
    end

    @testset "2D projection with different density and internal energy floors" begin
        equations = CompressibleEulerEquations2D(1.4)
        u = SVector(0.5, 1.0, -2.0, 0.1)
        lower_bounds = (1.0, 0.1)
        variables = (density, energy_internal)

        u_projected = Trixi.project_to_admissible_set(u, lower_bounds, variables, equations)
        arithmetic_tol = Trixi.euler_arithmetic_tol(lower_bounds[1], lower_bounds[2])

        @test u_projected[1] > lower_bounds[1]
        @test energy_internal(u_projected, equations) > lower_bounds[2] - arithmetic_tol
    end

    # this test failed without the rationalized approximation of
    # 0.5 * (rho - sqrt_discriminant_rho) in PositivityPreservingLimiterLiuZhang.
    @testset "2D projection near zero momentum boundary" begin
        equations = CompressibleEulerEquations2D(1.4)
        u = SVector(0.9376339560775117,
                    1.353902558446827e-8,
                    6.230911048510285e-9,
                    -3.779101287247884)
        lower_bounds = (1e-8, 2.5e-8)
        variables = (density, energy_internal)

        u_proj = Trixi.project_to_admissible_set(u, lower_bounds, variables, equations)
        arithmetic_tol = Trixi.euler_arithmetic_tol(lower_bounds[1], lower_bounds[2])

        @test u_proj[1]≈0.9376339560775118 rtol=1e-12
        @test u_proj[2]≈6.637197729990257e-9 rtol=1e-12
        @test u_proj[3]≈3.05456167498393e-9 rtol=1e-12
        @test u_proj[4]≈2.5000000028466725e-8 rtol=1e-12
        @test energy_internal(u_proj, equations) >= lower_bounds[2] - arithmetic_tol
    end
end
