using Printf, BenchmarkTools, Trixi

function run_benchmarks(benchmark_run; levels = 0:5, polydeg = 3)
    runtimes = zeros(length(levels))
    for (idx, initial_refinement_level) in enumerate(levels)
        result = benchmark_run(; initial_refinement_level, polydeg)
        display(result)
        runtimes[idx] = result |> median |> time # in nanoseconds
    end
    return (; levels, runtimes, polydeg)
end

function tabulate_benchmarks(args...; kwargs...)
    result = run_benchmarks(args...; kwargs...)
    println("#Elements | Runtime in seconds")
    for (level, runtime) in zip(result.levels, result.runtimes)
        @printf("%9d | %.2e\n", 4^level, 1.0e-9*runtime)
    end
    for (level, runtime) in zip(result.levels, result.runtimes)
        @printf("%.16e\n", 1.0e-9*runtime)
    end
end

function benchmark_euler(; initial_refinement_level = 1, polydeg = 3)
    γ = 1.4
    equations = CompressibleEulerEquations2D(γ)

    surface_flux = flux_ranocha
    volume_flux = flux_ranocha
    solver = DGSEM(polydeg, surface_flux, VolumeIntegralFluxDifferencing(volume_flux))

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

    semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_weak_blast_wave,
                                        solver;
                                        boundary_conditions = boundary_condition_periodic)

    t0 = 0.0
    u0 = compute_coefficients(t0, semi)
    du = similar(u0)

    @benchmark Trixi.rhs_hyperbolic!($du, $u0, $semi, $t0)
end

# versioninfo(verbose=true)
@show Threads.nthreads()
tabulate_benchmarks(benchmark_euler, levels = 0:8)
