@testsnippet Visualization begin
    using Plots
    # We use CairoMakie to avoid some CI-related issues with GLMakie. CairoMakie does not support
    # interactive visualization through `iplot`, but it can be used as a testing backend for Trixi's
    # Makie-based visualization.
    using CairoMakie

    EXAMPLES_DIR = examples_dir()
end

@testitem "Visualization: PlotData2D, PlotDataSeries, PlotMesh" setup=[Setup, Visualization] tags=[:misc_part1] begin

    # Run 2D tests with elixirs for different mesh and solver types
    test_examples_2d = Dict("TreeMesh" => ("tree_2d_dgsem",
                                           "elixir_euler_blast_wave_amr.jl"),
                            "TreeMesh (FDSBP)" => ("tree_2d_fdsbp",
                                                   "elixir_euler_convergence.jl"),
                            "StructuredMesh" => ("structured_2d_dgsem",
                                                 "elixir_euler_source_terms_waving_flag.jl"),
                            "UnstructuredMesh" => ("unstructured_2d_dgsem",
                                                   "elixir_euler_basic.jl"),
                            "P4estMesh" => ("p4est_2d_dgsem",
                                            "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"),
                            "DGMulti" => ("dgmulti_2d", "elixir_euler_weakform.jl"),
                            "DGMulti (FDSBP)" => ("dgmulti_2d",
                                                  "elixir_euler_cgsbp_periodic.jl"))

    @testset "PlotData2D, PlotDataSeries, PlotMesh with $mesh" for mesh in keys(test_examples_2d)
        # Run Trixi.jl
        directory, elixir = test_examples_2d[mesh]
        @test_trixi_include(joinpath(EXAMPLES_DIR, directory, elixir),
                            tspan=(0, 0.1))

        # Constructor tests
        if mesh == "TreeMesh" || mesh == "TreeMesh (FDSBP)"
            @test PlotData2D(sol) isa Trixi.PlotData2DCartesian
            @test PlotData2D(sol; nvisnodes = 0, grid_lines = false,
                             solution_variables = cons2cons) isa Trixi.PlotData2DCartesian
            if semi.solver isa DGSEM
                @test Trixi.PlotData2DTriangulated(sol) isa Trixi.PlotData2DTriangulated
            end
        else
            @test PlotData2D(sol) isa Trixi.PlotData2DTriangulated
            @test PlotData2D(sol; nvisnodes = 0, solution_variables = cons2cons) isa
                  Trixi.PlotData2DTriangulated
        end
        pd = PlotData2D(sol)

        # show
        @trixi_test_nowarn show(stdout, pd)
        println(stdout)

        # getindex
        @test pd["rho"] == Trixi.PlotDataSeries(pd, 1)
        @test pd["v1"] == Trixi.PlotDataSeries(pd, 2)
        @test pd["v2"] == Trixi.PlotDataSeries(pd, 3)
        @test pd["p"] == Trixi.PlotDataSeries(pd, 4)
        @test_throws KeyError pd["does not exist"]

        # convenience methods for mimicking a dictionary
        @test pd[begin] == Trixi.PlotDataSeries(pd, 1)
        @test pd[end] == Trixi.PlotDataSeries(pd, 4)
        @test length(pd) == 4
        @test size(pd) == (4,)
        @test keys(pd) == ("rho", "v1", "v2", "p")
        @test eltype(pd) <: Pair{String, <:Trixi.PlotDataSeries}
        @test [v for v in pd] == ["rho" => Trixi.PlotDataSeries(pd, 1),
            "v1" => Trixi.PlotDataSeries(pd, 2),
            "v2" => Trixi.PlotDataSeries(pd, 3),
            "p" => Trixi.PlotDataSeries(pd, 4)]

        # PlotDataSeries
        pds = pd["p"]
        @test pds.plot_data == pd
        @test pds.variable_id == 4
        @trixi_test_nowarn show(stdout, pds)
        println(stdout)

        # getmesh/PlotMesh
        @test getmesh(pd) == Trixi.PlotMesh(pd)
        @test getmesh(pd).plot_data == pd
        @trixi_test_nowarn show(stdout, getmesh(pd))
        println(stdout)

        @testset "2D plot recipes" begin
            pd = PlotData2D(sol)

            @trixi_test_nowarn Plots.plot(sol)
            @trixi_test_nowarn Plots.plot(pd)
            @trixi_test_nowarn Plots.plot(pd["p"])
            @trixi_test_nowarn Plots.plot(getmesh(pd))

            semi = sol.prob.p
            if mesh == "DGMulti"
                if sol.u[end] isa Trixi.VectorOfArray
                    u = parent(sol.u[end])
                else
                    u = sol.u[end]
                end
                scalar_data = StructArrays.component(u, 1)
            else
                # There are some issues with `PtrArray`s returned by default
                # by `Trixi.wrap_array`, see
                # https://github.com/trixi-framework/Trixi.jl/issues/2797
                u = Trixi.wrap_array_native(sol.u[end], semi)
                scalar_data = u[1, :, :, :]
            end

            if mesh != "DGMulti (FDSBP)"
                @trixi_test_nowarn Plots.plot(ScalarPlotData2D(scalar_data, semi))
                @trixi_test_nowarn Plots.plot(ScalarPlotData2D((u, equations) -> u[1],
                                                               sol.u[end], semi))
            end

            # test for consistency between the two ScalarPlotData2D constructions
            if mesh == "TreeMesh" || mesh == "TreeMesh (FDSBP)" || mesh == "DGMulti"
                spd_no_function = ScalarPlotData2D(scalar_data, semi)
                spd_function = ScalarPlotData2D((u, equations) -> u[1],
                                                sol.u[end], semi)
                @test typeof(spd_no_function) == typeof(spd_function)
                for property in propertynames(spd_function)
                    if property == :data
                        @test spd_no_function.data.data ≈ spd_function.data.data
                    elseif property == :variable_names
                        @test getproperty(spd_no_function, property) ==
                              getproperty(spd_function, property)
                    else
                        @test getproperty(spd_no_function, property) ≈
                              getproperty(spd_function, property)
                    end
                end
            end

            # test that non-upwinded FDSBP functionality is consistent with upwinded FDSBP solver behavior
            if mesh == "TreeMesh (FDSBP)"
                (; volume_integral, surface_integral) = semi.solver
                solver_fdsbp = FDSBP(semi.solver.basis.central; surface_integral,
                                     volume_integral)
                @test all(Trixi.reference_node_coordinates_2d(semi.solver) .≈
                          Trixi.reference_node_coordinates_2d(solver_fdsbp))
            end
        end

        @testset "1D plot from 2D solution" begin
            if mesh != "DGMulti" && mesh != "DGMulti (FDSBP)"
                @testset "Create 1D plot as slice" begin
                    @trixi_test_nowarn PlotData1D(sol, slice = :y, point = (0.5, 0.0)) isa
                                       PlotData1D
                    @trixi_test_nowarn PlotData1D(sol, slice = :x, point = (0.5, 0.0)) isa
                                       PlotData1D
                    pd1D = PlotData1D(sol, slice = :y, point = (0.5, 0.0))
                    @trixi_test_nowarn Plots.plot(pd1D)

                    @testset "Create 1D plot along curve" begin
                        curve = zeros(2, 10)
                        curve[1, :] = range(-1, 1, length = 10)
                        @trixi_test_nowarn PlotData1D(sol, curve = curve) isa PlotData1D
                        pd1D = PlotData1D(sol, curve = curve)
                        @trixi_test_nowarn Plots.plot(pd1D)
                    end
                end
            end
        end
    end
end
@testitem "Visualization: ScalarPlotData2D with DGMulti Quad elements" setup=[
    Setup,
    Visualization
] tags=[:misc_part1] begin

    # check that ScalarPlotData2D works for Quad elements
    @test_trixi_include(joinpath(EXAMPLES_DIR, "dgmulti_2d",
                                 "elixir_euler_weakform.jl"),
                        tspan=(0.0, 0.0),
                        element_type=Quad(),
                        cells_per_dimension=(2, 2))
    semi = sol.prob.p
    u = parent(sol.u[end])
    scalar_data = StructArrays.component(u, 1)
    @trixi_test_nowarn Plots.plot(ScalarPlotData2D(scalar_data, semi))
end
@testitem "Visualization: PlotData1D, PlotDataSeries, PlotMesh" setup=[Setup, Visualization] tags=[:misc_part1] begin

    # Run Trixi.jl
    @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem",
                                 "elixir_euler_blast_wave.jl"),
                        tspan=(0, 0.1))

    # Constructor
    @test PlotData1D(sol) isa PlotData1D
    pd = PlotData1D(sol)

    # show
    @trixi_test_nowarn show(stdout, pd)
    println(stdout)

    # getindex
    @test pd["rho"] == Trixi.PlotDataSeries(pd, 1)
    @test pd["v1"] == Trixi.PlotDataSeries(pd, 2)
    @test pd["p"] == Trixi.PlotDataSeries(pd, 3)
    @test_throws KeyError pd["does not exist"]

    # convenience methods for mimicking a dictionary
    @test pd[begin] == Trixi.PlotDataSeries(pd, 1)
    @test pd[end] == Trixi.PlotDataSeries(pd, 3)
    @test length(pd) == 3
    @test size(pd) == (3,)
    @test keys(pd) == ("rho", "v1", "p")
    @test eltype(pd) <: Pair{String, <:Trixi.PlotDataSeries}
    @test [v for v in pd] == ["rho" => Trixi.PlotDataSeries(pd, 1),
        "v1" => Trixi.PlotDataSeries(pd, 2),
        "p" => Trixi.PlotDataSeries(pd, 3)]

    # PlotDataSeries
    pds = pd["p"]
    @test pds.plot_data == pd
    @test pds.variable_id == 3
    @trixi_test_nowarn show(stdout, pds)
    println(stdout)

    # getmesh/PlotMesh
    @test getmesh(pd) == Trixi.PlotMesh(pd)
    @test getmesh(pd).plot_data == pd
    @trixi_test_nowarn show(stdout, getmesh(pd))
    println(stdout)

    # nvisnodes
    @test size(pd.data) == (512, 3)
    pd0 = PlotData1D(sol, nvisnodes = 0)
    @test size(pd0.data) == (256, 3)
    pd2 = PlotData1D(sol, nvisnodes = 2)
    @test size(pd2.data) == (128, 3)

    @testset "1D plot recipes" begin
        pd = PlotData1D(sol)

        @trixi_test_nowarn Plots.plot(sol)
        @trixi_test_nowarn Plots.plot(sol, reinterpolate = false)
        @trixi_test_nowarn Plots.plot(pd)
        @trixi_test_nowarn Plots.plot(pd["p"])
        @trixi_test_nowarn Plots.plot(getmesh(pd))
        initial_condition_t_end(x, equations) = initial_condition(x, last(tspan),
                                                                  equations)
        @trixi_test_nowarn Plots.plot(initial_condition_t_end, semi)
        @trixi_test_nowarn Plots.plot((x, equations) -> x, semi)
    end

    # Fake a PlotDataXD objects to test code for plotting multiple variables on at least two rows
    # with at least one plot remaining empty
    @testset "plotting multiple variables" begin
        x = collect(0.0:0.1:1.0)
        data1d = rand(5, 11)
        variable_names = string.('a':'e')
        mesh_vertices_x1d = [x[begin], x[end]]
        fake1d = PlotData1D(x, data1d, variable_names, mesh_vertices_x1d, 0)
        @trixi_test_nowarn Plots.plot(fake1d)

        y = x
        data2d = [rand(11, 11) for _ in 1:5]
        mesh_vertices_x2d = [0.0, 1.0, 1.0, 0.0]
        mesh_vertices_y2d = [0.0, 0.0, 1.0, 1.0]
        fake2d = Trixi.PlotData2DCartesian(x, y, data2d, variable_names,
                                           mesh_vertices_x2d, mesh_vertices_y2d, 0, 0)
        @trixi_test_nowarn Plots.plot(fake2d)
    end
end
@testitem "Visualization: FV testsets" setup=[Setup, Visualization] tags=[:misc_part1] begin
    @trixi_testset "BlockFV 1D Visualization" begin
        @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_blockfv",
                                     "elixir_advection_basic.jl"),
                            n_nodes=4,
                            initial_refinement_level=3)

        pd_fv_wrapped = @test_nowarn PlotData1D(sol)
        @test pd_fv_wrapped isa PlotData1D
    end

    @trixi_testset "DGSEM vs BlockFV 1D Visualization" begin
        @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem",
                                     "elixir_advection_finite_volume.jl"),
                            polydeg=0,
                            initial_refinement_level=5)

        pd_dgsem_polydeg0 = @test_nowarn PlotData1D(sol)

        @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_blockfv",
                                     "elixir_advection_basic.jl"),
                            n_nodes=4,
                            initial_refinement_level=3)

        pd_blockfv = @test_nowarn PlotData1D(sol)

        @test pd_blockfv.data ≈ pd_dgsem_polydeg0.data
    end

    @trixi_testset "Constant IC" begin
        @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_blockfv",
                                     "elixir_euler_source_term_nonperiodic.jl"),
                            n_nodes=4,
                            initial_refinement_level=3,
                            initial_condition=initial_condition_constant,
                            tspan=(0.0, 0.0))

        pd = @test_nowarn PlotData1D(sol)
        @test pd isa PlotData1D

        ref_cons = initial_condition_constant(SVector(0.0), 0.0,
                                              semi.equations)
        ref_prim = cons2prim(ref_cons, semi.equations)

        @test all(x -> isapprox(x, ref_prim[1]), pd.data[:, 1]) # rho
        @test all(x -> isapprox(x, ref_prim[2]), pd.data[:, 2]) # v1
        @test all(x -> isapprox(x, ref_prim[3]), pd.data[:, 3]) # p
    end
end
@testitem "Visualization: 1D plot from 2D solution" setup=[Setup, Visualization] tags=[:misc_part1] begin
    @trixi_testset "Create 1D plot along curve" begin
        using OrdinaryDiffEqSSPRK
        using Trixi

        @testset "$MeshType" for MeshType in (P4estMesh, T8codeMesh)
            equations = CompressibleEulerEquations2D(1.4)
            solver = DGSEM(polydeg = 3,
                           surface_flux = FluxLaxFriedrichs(max_abs_speed_naive))

            coordinates_min = (-1.0, -1.0)
            coordinates_max = (+1.0, +1.0)
            trees_per_dimension = (2, 2)
            mesh = MeshType(trees_per_dimension; polydeg = 1,
                            coordinates_min, coordinates_max,
                            initial_refinement_level = 0,
                            periodicity = true)

            semi = SemidiscretizationHyperbolic(mesh, equations,
                                                initial_condition_constant,
                                                solver;
                                                boundary_conditions = boundary_condition_periodic)
            ode = semidiscretize(semi, (0.0, 0.1))
            # SSPRK43 with optimized controller of Ranocha, Dalcin, Parsani,
            # and Ketcheson (2021)
            sol = solve(ode, SSPRK43();
                        dt = 0.01, ode_default_options()...)

            x = range(0, 1, length = 100)
            curve = vcat(x', x')
            # Ensure that PlotData1D is type stable
            pd = @inferred PlotData1D(sol; curve)

            # Compare with reference data
            u_node = SVector(1.0, 0.1, -0.2, 10.0)
            val = cons2prim(u_node, equations)
            for v in axes(pd.data, 2), i in axes(pd.data, 1)
                @test isapprox(pd.data[i, v], val[v])
            end
        end
    end

    @trixi_testset "PlotData1D gives correct results" begin
        using Trixi
        equations = CompressibleEulerEquations2D(1.4)
        solver = DGSEM(polydeg = 3,
                       surface_flux = FluxLaxFriedrichs(max_abs_speed_naive))
        coordinates_min = (-1.0, -1.0)
        coordinates_max = (+1.0, +1.0)
        initial_refinement_level = 3

        mesh_tree = TreeMesh(coordinates_min, coordinates_max;
                             initial_refinement_level,
                             periodicity = true)
        trees_per_dimension = (1, 1)
        mesh_p4est = P4estMesh(trees_per_dimension; polydeg = 1,
                               coordinates_min, coordinates_max,
                               initial_refinement_level,
                               periodicity = true)
        mesh_t8code = T8codeMesh(trees_per_dimension; polydeg = 1,
                                 coordinates_min, coordinates_max,
                                 initial_refinement_level,
                                 periodicity = true)
        cells_per_dimension = (2, 2) .^ initial_refinement_level
        mesh_structured = StructuredMesh(cells_per_dimension,
                                         coordinates_min, coordinates_max,
                                         periodicity = true)

        function initial_condition_taylor_green_vortex(x, t,
                                                       equations::CompressibleEulerEquations2D)
            A = 1.0 # magnitude of speed
            Ms = 0.1 # maximum Mach number

            rho = 1.0
            v1 = A * sin(x[1]) * cos(x[2])
            v2 = -A * cos(x[1]) * sin(x[2])
            p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms
            p = p +
                1.0 / 16.0 * A^2 * rho *
                (cos(2 * x[1]) + 2 * cos(2 * x[2]) +
                 2 * cos(2 * x[1]) + cos(2 * x[2]))

            return prim2cons(SVector(rho, v1, v2, p), equations)
        end
        ic = initial_condition_taylor_green_vortex

        ode_tree = let mesh = mesh_tree
            semi = SemidiscretizationHyperbolic(mesh, equations, ic, solver;
                                                boundary_conditions = boundary_condition_periodic)
            ode = semidiscretize(semi, (0.0, 0.1))
        end

        ode_p4est = let mesh = mesh_p4est
            semi = SemidiscretizationHyperbolic(mesh, equations, ic, solver;
                                                boundary_conditions = boundary_condition_periodic)
            ode = semidiscretize(semi, (0.0, 0.1))
        end

        ode_t8code = let mesh = mesh_t8code
            semi = SemidiscretizationHyperbolic(mesh, equations, ic, solver;
                                                boundary_conditions = boundary_condition_periodic)
            ode = semidiscretize(semi, (0.0, 0.1))
        end

        ode_structured = let mesh = mesh_structured
            semi = SemidiscretizationHyperbolic(mesh, equations, ic, solver;
                                                boundary_conditions = boundary_condition_periodic)
            ode = semidiscretize(semi, (0.0, 0.1))
        end

        x_curve = range(0, 1, length = 1000)
        curve = vcat(x_curve', x_curve')

        @testset "TreeMesh" begin
            pd_tree = @inferred PlotData1D(ode_tree.u0, ode_tree.p; curve)
            @test pd_tree.x ≈ range(0, 1, length = length(x_curve)) * sqrt(2)

            for i in eachindex(pd_tree.x)
                x = SVector(curve[1, i],
                            curve[2, i])
                prim = SVector(pd_tree.data[i, 1],
                               pd_tree.data[i, 2],
                               pd_tree.data[i, 3],
                               pd_tree.data[i, 4])
                u = initial_condition_taylor_green_vortex(x, 0.0, equations)
                @test isapprox(prim, cons2prim(u, equations), atol = 1.0e-3)
            end
        end

        @testset "P4estMesh" begin
            pd_p4est = @inferred PlotData1D(ode_p4est.u0, ode_p4est.p; curve)
            @test pd_p4est.x ≈ range(0, 1, length = length(x_curve)) * sqrt(2)

            for i in eachindex(pd_p4est.x)
                x = SVector(curve[1, i],
                            curve[2, i])
                prim = SVector(pd_p4est.data[i, 1],
                               pd_p4est.data[i, 2],
                               pd_p4est.data[i, 3],
                               pd_p4est.data[i, 4])
                u = initial_condition_taylor_green_vortex(x, 0.0, equations)
                @test isapprox(prim, cons2prim(u, equations), atol = 1.0e-3)
            end
        end

        @testset "T8codeMesh" begin
            pd_t8code = @inferred PlotData1D(ode_t8code.u0, ode_t8code.p; curve)
            @test pd_t8code.x ≈ range(0, 1, length = length(x_curve)) * sqrt(2)

            for i in eachindex(pd_t8code.x)
                x = SVector(curve[1, i],
                            curve[2, i])
                prim = SVector(pd_t8code.data[i, 1],
                               pd_t8code.data[i, 2],
                               pd_t8code.data[i, 3],
                               pd_t8code.data[i, 4])
                # Note that the T8codeMesh uses a different algorithm to
                # compute the 1D data than the other meshes above. This is
                # less accurate so that we need to use a larger tolerance.
                u = initial_condition_taylor_green_vortex(x, 0.0, equations)
                @test isapprox(prim, cons2prim(u, equations), atol = 5.0e-3)
            end
        end

        @testset "StructuredMesh" begin
            pd_structured = @inferred PlotData1D(ode_structured.u0, ode_structured.p;
                                                 curve)
            @test pd_structured.x ≈ range(0, 1, length = length(x_curve)) * sqrt(2)

            for i in eachindex(pd_structured.x)
                x = SVector(curve[1, i],
                            curve[2, i])
                prim = SVector(pd_structured.data[i, 1],
                               pd_structured.data[i, 2],
                               pd_structured.data[i, 3],
                               pd_structured.data[i, 4])
                # Note that the StructuredMesh uses a different algorithm to
                # compute the 1D data than the other meshes above. This is
                # less accurate so that we need to use a larger tolerance.
                u = initial_condition_taylor_green_vortex(x, 0.0, equations)
                @test isapprox(prim, cons2prim(u, equations), atol = 5.0e-3)
            end
        end
    end
end
@testitem "Visualization: PlotData2D Regression Tests" setup=[Setup, Visualization] tags=[:misc_part1] begin
    using Trixi
    equations = CompressibleEulerEquations2D(1.4)
    solver = DGSEM(polydeg = 3,
                   surface_flux = FluxLaxFriedrichs(max_abs_speed_naive))

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

    # Manually initialize meshes
    mesh_tree = TreeMesh(coordinates_min, coordinates_max;
                         initial_refinement_level,
                         periodicity = true)

    trees_per_dimension = (1, 1)
    mesh_p4est = P4estMesh(trees_per_dimension; polydeg = 3,
                           coordinates_min, coordinates_max,
                           initial_refinement_level,
                           periodicity = true)

    cells_per_dimension = (2, 2) .^ initial_refinement_level
    mesh_structured = StructuredMesh(cells_per_dimension,
                                     coordinates_min, coordinates_max,
                                     periodicity = true)

    function initial_condition_taylor_green_vortex(x, t,
                                                   equations::CompressibleEulerEquations2D)
        A = 1.0 # magnitude of speed
        Ms = 0.1 # maximum Mach number

        rho = 1.0
        v1 = A * sin(x[1]) * cos(x[2])
        v2 = -A * cos(x[1]) * sin(x[2])
        p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms
        p = p +
            1.0 / 16.0 * A^2 * rho *
            (cos(2 * x[1]) + 2 * cos(2 * x[2]) +
             2 * cos(2 * x[1]) + cos(2 * x[2]))

        return prim2cons(SVector(rho, v1, v2, p), equations)
    end

    @testset "Constant IC (Exact Checks)" begin
        ic = initial_condition_constant

        @testset "TreeMesh" begin
            semi_tree = SemidiscretizationHyperbolic(mesh_tree, equations, ic, solver;
                                                     boundary_conditions = boundary_condition_periodic)
            u_ode = compute_coefficients(0.0, semi_tree)
            pd = PlotData2D(u_ode, semi_tree, solution_variables = cons2prim)

            ref_cons = Trixi.initial_condition_constant(SVector(0.0, 0.0), 0.0,
                                                        semi_tree.equations)
            ref_prim = cons2prim(ref_cons, semi_tree.equations)

            @test all(x -> isapprox(x, ref_prim[1]), pd.data[1]) # rho
            @test all(x -> isapprox(x, ref_prim[2]), pd.data[2]) # v1
            @test all(x -> isapprox(x, ref_prim[3]), pd.data[3]) # v2
            @test all(x -> isapprox(x, ref_prim[4]), pd.data[4]) # p
        end

        @testset "StructuredMesh" begin
            semi_struct = SemidiscretizationHyperbolic(mesh_structured, equations, ic,
                                                       solver;
                                                       boundary_conditions = boundary_condition_periodic)
            u_ode = compute_coefficients(0.0, semi_struct)
            pd = PlotData2D(u_ode, semi_struct, solution_variables = cons2prim)

            ref_cons = Trixi.initial_condition_constant(SVector(0.0, 0.0), 0.0,
                                                        semi_struct.equations)
            ref_prim = cons2prim(ref_cons, semi_struct.equations)

            @test all(val -> isapprox(val[1], ref_prim[1]), pd.data) # rho
            @test all(val -> isapprox(val[2], ref_prim[2]), pd.data) # v1
            @test all(val -> isapprox(val[3], ref_prim[3]), pd.data) # v2
            @test all(val -> isapprox(val[4], ref_prim[4]), pd.data) # p
        end

        @testset "P4estMesh" begin
            semi_p4est = SemidiscretizationHyperbolic(mesh_p4est, equations, ic, solver;
                                                      boundary_conditions = boundary_condition_periodic)
            u_ode = compute_coefficients(0.0, semi_p4est)
            pd = PlotData2D(u_ode, semi_p4est, solution_variables = cons2prim)

            ref_cons = Trixi.initial_condition_constant(SVector(0.0, 0.0), 0.0,
                                                        semi_p4est.equations)
            ref_prim = cons2prim(ref_cons, semi_p4est.equations)

            @test all(val -> isapprox(val[1], ref_prim[1]), pd.data) # rho
            @test all(val -> isapprox(val[2], ref_prim[2]), pd.data) # v1
            @test all(val -> isapprox(val[3], ref_prim[3]), pd.data) # v2
            @test all(val -> isapprox(val[4], ref_prim[4]), pd.data) # p
        end
    end

    @testset "Non-Constant IC (Taylor-Green Vortex)" begin
        ic = initial_condition_taylor_green_vortex

        @testset "TreeMesh" begin
            semi_tree = SemidiscretizationHyperbolic(mesh_tree, equations, ic, solver;
                                                     boundary_conditions = boundary_condition_periodic)
            u_ode = compute_coefficients(0.0, semi_tree)
            pd = PlotData2D(u_ode, semi_tree, solution_variables = cons2prim)

            max_error = 0.0
            for (j, y) in enumerate(pd.y), (i, x) in enumerate(pd.x)
                u_exact = ic(SVector(x, y), 0.0, semi_tree.equations)
                prim_exact = cons2prim(u_exact, semi_tree.equations)
                prim_interp = SVector(pd.data[1][i, j], pd.data[2][i, j],
                                      pd.data[3][i, j], pd.data[4][i, j])

                current_error = maximum(abs.(prim_interp - prim_exact))
                max_error = max(max_error, current_error)
            end
            # Note that PlotData2D for TreeMesh uses a different algorithm that interpolates
            # the solution onto a uniform Cartesian grid. This is less accurate than the
            # exact nodal evaluations above, so we need to use a larger tolerance.
            @test max_error < 1.05
        end

        @testset "StructuredMesh" begin
            semi_struct = SemidiscretizationHyperbolic(mesh_structured, equations, ic,
                                                       solver;
                                                       boundary_conditions = boundary_condition_periodic)
            u_ode = compute_coefficients(0.0, semi_struct)
            pd = PlotData2D(u_ode, semi_struct, solution_variables = cons2prim)

            max_error = 0.0
            for i in eachindex(pd.x)
                x = pd.x[i]
                y = pd.y[i]
                u_exact = ic(SVector(x, y), 0.0, semi_struct.equations)
                prim_exact = cons2prim(u_exact, semi_struct.equations)

                current_error = maximum(abs.(pd.data[i] - prim_exact))
                max_error = max(max_error, current_error)
            end
            @test max_error < 1.0e-5
        end

        @testset "P4estMesh" begin
            semi_p4est = SemidiscretizationHyperbolic(mesh_p4est, equations, ic, solver;
                                                      boundary_conditions = boundary_condition_periodic)
            u_ode = compute_coefficients(0.0, semi_p4est)
            pd = PlotData2D(u_ode, semi_p4est, solution_variables = cons2prim)

            max_error = 0.0
            for i in eachindex(pd.x)
                x = pd.x[i]
                y = pd.y[i]
                u_exact = ic(SVector(x, y), 0.0, semi_p4est.equations)
                prim_exact = cons2prim(u_exact, semi_p4est.equations)

                current_error = maximum(abs.(pd.data[i] - prim_exact))
                max_error = max(max_error, current_error)
            end
            @test max_error < 1.0e-5
        end
    end
end
@testitem "Visualization: PlotData1D (DGMulti)" setup=[Setup, Visualization] tags=[:misc_part1] begin

    # Test two different approximation types since these use different memory layouts:
    # - structure of arrays for `Polynomial()`
    # - array of structures for `SBP()`
    @test_trixi_include(joinpath(EXAMPLES_DIR, "dgmulti_1d",
                                 "elixir_euler_flux_diff.jl"),
                        tspan=(0.0, 0.0),
                        approximation_type=Polynomial())
    @test PlotData1D(sol) isa PlotData1D
    initial_condition_t_end(x, equations) = initial_condition(x, last(tspan), equations)
    @trixi_test_nowarn Plots.plot(initial_condition_t_end, semi)
    @trixi_test_nowarn Plots.plot((x, equations) -> x, semi)

    @test_trixi_include(joinpath(EXAMPLES_DIR, "dgmulti_1d",
                                 "elixir_euler_flux_diff.jl"),
                        tspan=(0.0, 0.0),
                        approximation_type=SBP())
    @test PlotData1D(sol) isa PlotData1D
    @trixi_test_nowarn Plots.plot(initial_condition_t_end, semi)
    @trixi_test_nowarn Plots.plot((x, equations) -> x, semi)
end
@testitem "Visualization: PlotData2D (DGMulti Tri SBP)" setup=[Setup, Visualization] tags=[:misc_part1] begin

    # Regression test for plotting with SBP on triangular elements (reference triangulation of rstp).
    @test_trixi_include(joinpath(EXAMPLES_DIR, "dgmulti_2d",
                                 "elixir_euler_weakform.jl"),
                        cells_per_dimension=(4, 4),
                        approximation_type=SBP(),
                        surface_integral=SurfaceIntegralWeakForm(FluxHLL(min_max_speed_naive)),
                        tspan=(0.0, 0.0))

    pd = PlotData2D(sol)
    @test pd isa Trixi.PlotData2DTriangulated
    @test size(pd.t, 1) > 0

    @trixi_test_nowarn Plots.plot(pd)
    @trixi_test_nowarn Plots.plot(pd["rho"])
end
@testitem "Visualization: 1D plot recipes (StructuredMesh)" setup=[Setup, Visualization] tags=[:misc_part1] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "structured_1d_dgsem",
                                 "elixir_euler_source_terms.jl"),
                        tspan=(0.0, 0.0))

    pd = PlotData1D(sol)
    initial_condition_t_end(x, equations) = initial_condition(x, last(tspan), equations)
    @trixi_test_nowarn Plots.plot(sol)
    @trixi_test_nowarn Plots.plot(pd)
    @trixi_test_nowarn Plots.plot(pd["p"])
    @trixi_test_nowarn Plots.plot(sol.u[end], semi)
    @trixi_test_nowarn Plots.plot(initial_condition_t_end, semi)
    @trixi_test_nowarn Plots.plot((x, equations) -> x, semi)
end
@testitem "Visualization: plot time series" setup=[Setup, Visualization] tags=[:misc_part1] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
                                 "elixir_acoustics_gaussian_source.jl"),
                        tspan=(0, 0.05))

    @trixi_test_nowarn Plots.plot(time_series, 1)
    @test PlotData1D(time_series, 1) isa PlotData1D
end
@testitem "Visualization: adapt_to_mesh_level" setup=[Setup, Visualization] tags=[:misc_part1] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
                                 "elixir_advection_basic.jl"),
                        analysis_callback=Trixi.TrivialCallback())
    @test adapt_to_mesh_level(sol, 5) isa Tuple

    u_ode_level5, semi_level5 = adapt_to_mesh_level(sol, 5)
    u_ode_level4, semi_level4 = adapt_to_mesh_level(u_ode_level5, semi_level5, 4)
    @test isapprox(sol.u[end], u_ode_level4, atol = 1e-13)

    @test adapt_to_mesh_level!(sol, 5) isa Tuple
    @test isapprox(sol.u[end], u_ode_level5, atol = 1e-13)
end
@testitem "Visualization: plot 3D" setup=[Setup, Visualization] tags=[:misc_part1] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_3d_dgsem",
                                 "elixir_advection_basic.jl"),
                        analysis_callback=Trixi.TrivialCallback(),
                        initial_refinement_level=1)
    @test PlotData2D(sol) isa Trixi.PlotData2DCartesian
    @test PlotData2D(sol, slice = :yz) isa Trixi.PlotData2DCartesian
    @test PlotData2D(sol, slice = :xz) isa Trixi.PlotData2DCartesian

    @testset "1D plot from 3D solution and Tree-mesh" begin
        @testset "Create 1D plot as slice" begin
            @trixi_test_nowarn PlotData1D(sol) isa PlotData1D
            pd1D = PlotData1D(sol)
            @trixi_test_nowarn Plots.plot(pd1D)
            @trixi_test_nowarn PlotData1D(sol, slice = :y, point = (0.5, 0.3, 0.1)) isa
                               PlotData1D
            @trixi_test_nowarn PlotData1D(sol, slice = :z, point = (0.1, 0.3, 0.3)) isa
                               PlotData1D
        end

        @testset "Create 1D plot along curve" begin
            curve = zeros(3, 10)
            curve[1, :] = range(-1.0, -0.5, length = 10)
            @trixi_test_nowarn PlotData1D(sol, curve = curve) isa PlotData1D
            pd1D = PlotData1D(sol, curve = curve)
            @trixi_test_nowarn Plots.plot(pd1D)
        end
    end

    @test_trixi_include(joinpath(EXAMPLES_DIR, "structured_3d_dgsem",
                                 "elixir_advection_basic.jl"))

    @testset "1D plot from 3D solution and general mesh" begin
        @testset "Create 1D plot as slice" begin
            @trixi_test_nowarn PlotData1D(sol) isa PlotData1D
            pd1D = PlotData1D(sol)
            @trixi_test_nowarn Plots.plot(pd1D)
            @trixi_test_nowarn PlotData1D(sol, slice = :y, point = (0.5, 0.3, 0.1)) isa
                               PlotData1D
            @trixi_test_nowarn PlotData1D(sol, slice = :z, point = (0.1, 0.3, 0.3)) isa
                               PlotData1D
        end

        @testset "Create 1D plot along curve" begin
            curve = zeros(3, 10)
            curve[1, :] = range(-1.0, 1.0, length = 10)
            @trixi_test_nowarn PlotData1D(sol, curve = curve) isa PlotData1D
            pd1D = PlotData1D(sol, curve = curve)
            @trixi_test_nowarn Plots.plot(pd1D)
        end
    end

    @timed_testset "1D plot from 3D solution on P4estMesh and T8codeMesh" begin
        @trixi_testset "Create 1D plot along curve" begin
            using OrdinaryDiffEqSSPRK
            using Trixi
            @testset "$MeshType" for MeshType in (P4estMesh, T8codeMesh)
                equations = CompressibleEulerEquations3D(1.4)
                solver = DGSEM(polydeg = 3,
                               surface_flux = FluxLaxFriedrichs(max_abs_speed_naive))

                coordinates_min = (-1.0, -1.0, -1.0)
                coordinates_max = (+1.0, +1.0, +1.0)
                trees_per_dimension = (2, 2, 2)
                mesh = MeshType(trees_per_dimension; polydeg = 1,
                                coordinates_min, coordinates_max,
                                initial_refinement_level = 0,
                                periodicity = true)

                semi = SemidiscretizationHyperbolic(mesh, equations,
                                                    initial_condition_constant,
                                                    solver;
                                                    boundary_conditions = boundary_condition_periodic)
                ode = semidiscretize(semi, (0.0, 0.1))
                # SSPRK43 with optimized controller of Ranocha, Dalcin, Parsani,
                # and Ketcheson (2021)
                sol = solve(ode, SSPRK43();
                            dt = 0.01, ode_default_options()...)

                x = range(0, 1, length = 100)
                curve = vcat(x', x', x')
                # Ensure that PlotData1D is type stable
                pd = @inferred PlotData1D(sol; curve)

                # Compare with reference data
                u_node = SVector(1.0, 0.1, -0.2, 0.7, 10.0)
                val = cons2prim(u_node, equations)
                for v in axes(pd.data, 2), i in axes(pd.data, 1)
                    @test isapprox(pd.data[i, v], val[v])
                end
            end
        end
    end

    @trixi_testset "PlotData1D gives correct results" begin
        using Trixi
        equations = CompressibleEulerEquations3D(1.4)
        solver = DGSEM(polydeg = 3,
                       surface_flux = FluxLaxFriedrichs(max_abs_speed_naive))
        coordinates_min = (-1.0, -1.0, -1.0)
        coordinates_max = (+1.0, +1.0, +1.0)
        initial_refinement_level = 3

        mesh_tree = TreeMesh(coordinates_min, coordinates_max;
                             initial_refinement_level,
                             periodicity = true)
        trees_per_dimension = (1, 1, 1)
        mesh_p4est = P4estMesh(trees_per_dimension; polydeg = 1,
                               coordinates_min, coordinates_max,
                               initial_refinement_level,
                               periodicity = true)
        mesh_t8code = T8codeMesh(trees_per_dimension; polydeg = 1,
                                 coordinates_min, coordinates_max,
                                 initial_refinement_level,
                                 periodicity = true)
        cells_per_dimension = (2, 2, 2) .^ initial_refinement_level
        mesh_structured = StructuredMesh(cells_per_dimension,
                                         coordinates_min, coordinates_max,
                                         periodicity = true)

        function initial_condition_taylor_green_vortex(x, t,
                                                       equations::CompressibleEulerEquations3D)
            A = 1.0 # magnitude of speed
            Ms = 0.1 # maximum Mach number

            rho = 1.0
            v1 = A * sin(x[1]) * cos(x[2]) * cos(x[3])
            v2 = -A * cos(x[1]) * sin(x[2]) * cos(x[3])
            v3 = 0.0
            p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms
            p = p +
                1.0 / 16.0 * A^2 * rho *
                (cos(2 * x[1]) * cos(2 * x[3]) + 2 * cos(2 * x[2]) +
                 2 * cos(2 * x[1]) +
                 cos(2 * x[2]) * cos(2 * x[3]))

            return prim2cons(SVector(rho, v1, v2, v3, p), equations)
        end
        ic = initial_condition_taylor_green_vortex

        ode_tree = let mesh = mesh_tree
            semi = SemidiscretizationHyperbolic(mesh, equations, ic, solver;
                                                boundary_conditions = boundary_condition_periodic)
            ode = semidiscretize(semi, (0.0, 0.1))
        end

        ode_p4est = let mesh = mesh_p4est
            semi = SemidiscretizationHyperbolic(mesh, equations, ic, solver;
                                                boundary_conditions = boundary_condition_periodic)
            ode = semidiscretize(semi, (0.0, 0.1))
        end

        ode_t8code = let mesh = mesh_t8code
            semi = SemidiscretizationHyperbolic(mesh, equations, ic, solver;
                                                boundary_conditions = boundary_condition_periodic)
            ode = semidiscretize(semi, (0.0, 0.1))
        end

        ode_structured = let mesh = mesh_structured
            semi = SemidiscretizationHyperbolic(mesh, equations, ic, solver;
                                                boundary_conditions = boundary_condition_periodic)
            ode = semidiscretize(semi, (0.0, 0.1))
        end

        x_curve = range(0, 1, length = 1000)
        curve = vcat(x_curve', x_curve', x_curve')

        @testset "TreeMesh" begin
            pd_tree = @inferred PlotData1D(ode_tree.u0, ode_tree.p; curve)
            @test pd_tree.x ≈ range(0, 1, length = length(x_curve)) * sqrt(3)

            for i in eachindex(pd_tree.x)
                x = SVector(curve[1, i],
                            curve[2, i],
                            curve[3, i])
                prim = SVector(pd_tree.data[i, 1],
                               pd_tree.data[i, 2],
                               pd_tree.data[i, 3],
                               pd_tree.data[i, 4],
                               pd_tree.data[i, 5])
                u = initial_condition_taylor_green_vortex(x, 0.0, equations)
                @test isapprox(prim, cons2prim(u, equations), atol = 1.0e-3)
            end
        end

        @testset "P4estMesh" begin
            pd_p4est = @inferred PlotData1D(ode_p4est.u0, ode_p4est.p; curve)
            @test pd_p4est.x ≈ range(0, 1, length = length(x_curve)) * sqrt(3)

            for i in eachindex(pd_p4est.x)
                x = SVector(curve[1, i],
                            curve[2, i],
                            curve[3, i])
                prim = SVector(pd_p4est.data[i, 1],
                               pd_p4est.data[i, 2],
                               pd_p4est.data[i, 3],
                               pd_p4est.data[i, 4],
                               pd_p4est.data[i, 5])
                u = initial_condition_taylor_green_vortex(x, 0.0, equations)
                @test isapprox(prim, cons2prim(u, equations), atol = 1.0e-3)
            end
        end

        @testset "T8codeMesh" begin
            pd_t8code = @inferred PlotData1D(ode_t8code.u0, ode_t8code.p; curve)
            @test pd_t8code.x ≈ range(0, 1, length = length(x_curve)) * sqrt(3)

            for i in eachindex(pd_t8code.x)
                x = SVector(curve[1, i],
                            curve[2, i],
                            curve[3, i])
                prim = SVector(pd_t8code.data[i, 1],
                               pd_t8code.data[i, 2],
                               pd_t8code.data[i, 3],
                               pd_t8code.data[i, 4],
                               pd_t8code.data[i, 5])
                u = initial_condition_taylor_green_vortex(x, 0.0, equations)
                @test isapprox(prim, cons2prim(u, equations), atol = 1.0e-3)
            end
        end

        @testset "StructuredMesh" begin
            pd_structured = @inferred PlotData1D(ode_structured.u0, ode_structured.p;
                                                 curve)
            @test pd_structured.x ≈ range(0, 1, length = length(x_curve)) * sqrt(3)

            for i in eachindex(pd_structured.x)
                x = SVector(curve[1, i],
                            curve[2, i],
                            curve[3, i])
                prim = SVector(pd_structured.data[i, 1],
                               pd_structured.data[i, 2],
                               pd_structured.data[i, 3],
                               pd_structured.data[i, 4],
                               pd_structured.data[i, 5])
                u = initial_condition_taylor_green_vortex(x, 0.0, equations)
                # Note that the StructuredMesh uses a different algorithm to
                # compute the 1D data than the other meshes. This is less accurate
                # so that we need to use a larger tolerance.
                @test isapprox(prim, cons2prim(u, equations), atol = 5.0e-3)
            end
        end
    end
end
@testitem "Visualization: plotting TimeIntegratorSolution" setup=[Setup, Visualization] tags=[:misc_part1] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
                                 "elixir_hypdiff_lax_friedrichs.jl"),
                        maxiters=1, analysis_callback=Trixi.TrivialCallback(),
                        initial_refinement_level=1)
    @trixi_test_nowarn Plots.plot(sol)
end
@testitem "Visualization: VisualizationCallback" setup=[Setup, Visualization] tags=[:misc_part1] begin
    outdir = "out"

    # To make CI tests work, disable showing a plot window with the GR backend of the Plots package
    # Xref: https://github.com/jheinen/GR.jl/issues/278
    # Xref: https://github.com/JuliaPlots/Plots.jl/blob/8cc6d9d48755ba452a2835f9b89d3880e9945377/test/runtests.jl#L103
    if !isinteractive()
        restore = get(ENV, "GKSwstype", nothing)
        ENV["GKSwstype"] = "100"
    end

    @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
                                 "elixir_advection_amr_visualization.jl"),
                        visualization=VisualizationCallback(semi;
                                                            interval = 20,
                                                            clims = (0, 1),
                                                            plot_creator = Trixi.save_plot),
                        tspan=(0.0, 3.0))

    @testset "elixir_advection_amr_visualization.jl with save_plot" begin
        @test isfile(joinpath(outdir, "solution_000000000.png"))
        @test isfile(joinpath(outdir, "solution_000000020.png"))
        @test isfile(joinpath(outdir, "solution_000000022.png"))
    end

    @testset "show" begin
        @trixi_test_nowarn show(stdout, visualization)
        println(stdout)

        @trixi_test_nowarn show(stdout, "text/plain", visualization)
        println(stdout)
    end

    # Restore GKSwstype to previous value (if it was set)
    if !isinteractive()
        if isnothing(restore)
            delete!(ENV, "GKSwstype")
        else
            ENV["GKSwstype"] = restore
        end
    end
end
@testitem "Visualization: Makie visualization tests for 1D" setup=[Setup, Visualization] tags=[:misc_part1] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem",
                                 "elixir_advection_basic.jl"))
    pd = PlotData1D(sol)

    # convert_arguments enables lines(pd["scalar"])
    @trixi_test_nowarn lines(pd["scalar"])

    # plottype for 1D PlotDataSeries is Lines
    @test Makie.plottype(pd["scalar"]) == Makie.Lines

    # Makie.plot(pds) gives title and xlabel as for Plots.jl recipes
    @trixi_test_nowarn Makie.plot(pd["scalar"])

    # plot_mesh kwarg triggers vlines!
    @trixi_test_nowarn Makie.plot(pd["scalar"], plot_mesh = true)

    # kwargs are forwarded to lines!
    @trixi_test_nowarn Makie.plot(pd["scalar"], color = :red, linewidth = 2)

    # Makie.plot(pd) gives layout for all variables
    fa = Makie.plot(pd)

    # plot_mesh kwarg triggers vlines! for mesh vertices
    @trixi_test_nowarn Makie.plot(pd, plot_mesh = true)
    fig, axes = fa
    @trixi_test_nowarn Base.show(fa) === nothing

    # Makie.plot(sol) for 1D solutions
    @trixi_test_nowarn Makie.plot(sol)

    # PlotMesh overlay
    Makie.plot(pd["scalar"])
    @trixi_test_nowarn Makie.plot!(Trixi.PlotMesh(pd))

    # kwargs are forwarded to vlines! in PlotMesh
    Makie.plot(pd["scalar"])
    @trixi_test_nowarn Makie.plot!(Trixi.PlotMesh(pd), color = :black,
                                   linestyle = :dash)
end
@testitem "Visualization: Makie visualization tests for TreeMesh2D" setup=[
    Setup,
    Visualization
] tags=[:misc_part1] begin
    using CairoMakie
    @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
                                 "elixir_advection_basic.jl"))
    pd = @inferred PlotData2D(sol)
    @test pd isa Trixi.PlotData2DCartesian

    # plottype for 2D PlotDataSeries is Heatmap
    @test Makie.plottype(pd["scalar"]) == Makie.Heatmap

    # convert_arguments enables Makie.heatmap(pd["scalar"])
    @trixi_test_nowarn Makie.heatmap(pd["scalar"])

    # Makie.plot(pds) gives title, xlabel, ylabel and colorbar
    @trixi_test_nowarn Makie.plot(pd["scalar"])

    # kwargs are forwarded to heatmap!
    @trixi_test_nowarn Makie.plot(pd["scalar"], colormap = :heat)

    # Makie.plot(pd) gives layout for all variables
    fa = @trixi_test_nowarn Makie.plot(pd)
    @trixi_test_nowarn Makie.plot(pd, plot_mesh = true)
    fig, axes = fa
    @trixi_test_nowarn Base.show(fa) === nothing

    # Makie.plot(sol) for 2D TreeMesh solutions
    @trixi_test_nowarn Makie.plot(sol)

    # PlotMesh overlay
    Makie.plot(pd["scalar"])
    @trixi_test_nowarn Makie.plot!(Trixi.PlotMesh(pd))

    # kwargs are forwarded to lines! in PlotMesh
    Makie.plot(pd["scalar"])
    @trixi_test_nowarn Makie.plot!(Trixi.PlotMesh(pd), color = :black,
                                   linestyle = :dash)
end
@testitem "Visualization: Makie visualization tests for UnstructuredMesh2D" setup=[
    Setup,
    Visualization
] tags=[:misc_part1] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "unstructured_2d_dgsem",
                                 "elixir_euler_wall_bc.jl"))

    # test interactive surface plot
    @trixi_test_nowarn Trixi.iplot(sol)

    # also test when using PlotData2D object
    @test PlotData2D(sol) isa Trixi.PlotData2DTriangulated
    @trixi_test_nowarn Makie.plot(PlotData2D(sol))

    # test interactive ScalarPlotData2D plotting
    semi = sol.prob.p
    x = view(semi.cache.elements.node_coordinates, 1, :, :, :) # extracts the node x coordinates
    y = view(semi.cache.elements.node_coordinates, 2, :, :, :) # extracts the node x coordinates
    @trixi_test_nowarn iplot(ScalarPlotData2D(x .+ y, semi), plot_mesh = true)

    # test heatmap plot
    @trixi_test_nowarn Makie.plot(sol, plot_mesh = true)

    # single-variable plot with axis and colorbar (works for all PlotData2DTriangulated meshes)
    pd = @inferred PlotData2D(sol)
    @trixi_test_nowarn Makie.plot(pd["rho"])
    @trixi_test_nowarn Makie.plot(pd["rho"], colormap = :blues)
    # plot_mesh = true
    @trixi_test_nowarn Makie.plot(pd["rho"], plot_mesh = true)

    # explicit PlotMesh overlay (works for all PlotData2DTriangulated meshes)
    @trixi_test_nowarn Makie.plot(pd["rho"])
    @trixi_test_nowarn Makie.plot!(getmesh(pd))
    @trixi_test_nowarn Makie.plot(pd["rho"])
    @trixi_test_nowarn Makie.plot!(getmesh(pd), color = :black,
                                   linestyle = :dash)

    # test unpacking/iteration for FigureAndAxes
    fa = Makie.plot(sol)
    fig, axes = fa
    @trixi_test_nowarn Base.show(fa) === nothing
    @trixi_test_nowarn typeof(fig) <: Makie.Figure
    @trixi_test_nowarn typeof(axes) <: AbstractArray{<:Makie.Axis}

    # test plotting of constant solutions with Makie
    # related issue: https://github.com/MakieOrg/Makie.jl/issues/931
    for i in eachindex(sol.u)
        fill!(sol.u[i], one(eltype(sol.u[i])))
    end
    @trixi_test_nowarn Trixi.iplot(sol)
end
@testitem "Visualization: Makie iplot for DGMulti with VectorOfArray solution" setup=[
    Setup,
    Visualization
] tags=[:misc_part1] begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "dgmulti_2d",
                                 "elixir_euler_curved.jl"))

    @trixi_test_nowarn Trixi.iplot(sol)
end
