Add benchmarks

Author: efaulhaber

benchmarks/benchmarks.jl

@@ -0,0 +1,4 @@
+include("count_neighbors.jl")
+include("n_body.jl")
+
+include("plot.jl")

benchmarks/count_neighbors.jl

@@ -0,0 +1,29 @@
+using PointNeighbors
+using BenchmarkTools
+
+"""
+    benchmark_count_neighbors(neighborhood_search, coordinates; parallel = true)
+
+A very cheap and simple neighborhood search benchmark, only counting the neighbors of each
+particle. For each particle-neighbor pair, only an array entry is incremented.
+
+Due to the minimal computational cost, differences between neighborhood search
+implementations are highlighted. On the other hand, this is the least realistic benchmark.
+
+For a computationally heavier benchmark, see [`benchmark_n_body`](@ref).
+"""
+function benchmark_count_neighbors(neighborhood_search, coordinates; parallel = true)
+    n_neighbors = zeros(Int, size(coordinates, 2))
+
+    function count_neighbors!(n_neighbors, coordinates, neighborhood_search, parallel)
+        n_neighbors .= 0
+
+        for_particle_neighbor(coordinates, coordinates, neighborhood_search,
+                              parallel = parallel) do i, _, _, _
+            n_neighbors[i] += 1
+        end
+    end
+
+    return @belapsed $count_neighbors!($n_neighbors, $coordinates,
+                                       $neighborhood_search, $parallel)
+end

benchmarks/n_body.jl

@@ -0,0 +1,41 @@
+using PointNeighbors
+using BenchmarkTools
+
+"""
+    benchmark_n_body(neighborhood_search, coordinates; parallel = true)
+
+A simple neighborhood search benchmark, computing the right-hand side of an n-body
+simulation with a cutoff (corresponding to the search radius of `neighborhood_search`).
+
+This is a more realistic benchmark for particle-based simulations than
+[`benchmark_count_neighbors`](@ref).
+However, due to the higher computational cost, differences between neighborhood search
+implementations are less pronounced.
+"""
+function benchmark_n_body(neighborhood_search, coordinates; parallel = true)
+    mass = 1e10 * (rand(size(coordinates, 2)) .+ 1)
+    G = 6.6743e-11
+
+    dv = similar(coordinates)
+
+    function compute_acceleration!(dv, coordinates, mass, G, neighborhood_search, parallel)
+        dv .= 0.0
+
+        for_particle_neighbor(coordinates, coordinates, neighborhood_search,
+                              parallel = parallel) do i, j, pos_diff, distance
+            # Only consider particles with a distance > 0
+            distance < sqrt(eps()) && return
+
+            dv_ = -G * mass[j] * pos_diff / distance^3
+
+            for dim in axes(dv, 1)
+                @inbounds dv[dim, i] += dv_[dim]
+            end
+        end
+
+        return dv
+    end
+
+    return @belapsed $compute_acceleration!($dv, $coordinates, $mass, $G,
+                                            $neighborhood_search, $parallel)
+end

benchmarks/plot.jl

@@ -0,0 +1,78 @@
+using Plots
+using BenchmarkTools
+
+# Generate a rectangular point cloud
+include("../test/point_cloud.jl")
+
+"""
+    plot_benchmarks(benchmark, n_points_per_dimension, iterations;
+                    seed = 1, perturbation_factor_position = 1.0,
+                    parallel = true, title = "")
+
+Run a benchmark for with several neighborhood searches multiple times for increasing numbers
+of points and plot the results.
+
+# Arguments
+- `benchmark`:              The benchmark function. See [`benchmark_count_neigbors`](@ref)
+                            and [`benchmark_n_body`](@ref).
+- `n_points_per_dimension`: Initial resolution as tuple. The product is the initial number
+                            of points. For example, use `(100, 100)` for a 2D benchmark or
+                            `(10, 10, 10)` for a 3D benchmark.
+- `iterations`:             Number of refinement iterations
+
+# Keywords
+- `parallel = true`:        Loop over all points in parallel
+- `title = ""`:             Title of the plot
+- `seed = 1`:               Seed to perturb the point positions. Different seeds yield
+                            slightly different point positions.
+- `perturbation_factor_position = 1.0`: Perturb point positions by this factor. A factor of
+                                        `1.0` corresponds to points being moved by
+                                        a maximum distance of `0.5` along each axis.
+"""
+function plot_benchmarks(benchmark, n_points_per_dimension, iterations;
+                         parallel = true, title = "",
+                         seed = 1, perturbation_factor_position = 1.0)
+    neighborhood_searches_names = ["TrivialNeighborhoodSearch";;
+                                   "GridNeighborhoodSearch"]
+
+    # Multiply number of points in each iteration (roughly) by this factor
+    scaling_factor = 4
+    per_dimension_factor = scaling_factor^(1 / length(n_points_per_dimension))
+    sizes = [round.(Int, n_points_per_dimension .* per_dimension_factor^(iter - 1))
+             for iter in 1:iterations]
+
+    n_particles = prod.(sizes)
+    times = zeros(iterations, length(neighborhood_searches_names))
+
+    for iter in 1:iterations
+        coordinates = point_cloud(sizes[iter], seed = seed,
+                                  perturbation_factor_position = perturbation_factor_position)
+
+        search_radius = 3.0
+
+        neighborhood_searches = [
+            TrivialNeighborhoodSearch{size(coordinates, 1)}(search_radius,
+                                                            axes(coordinates, 2)),
+            GridNeighborhoodSearch{size(coordinates, 1)}(search_radius,
+                                                         size(coordinates, 2)),
+        ]
+
+        for i in eachindex(neighborhood_searches)
+            neighborhood_search = neighborhood_searches[i]
+            initialize!(neighborhood_search, coordinates, coordinates)
+
+            time = benchmark(neighborhood_search, coordinates, parallel = parallel)
+            times[iter, i] = time
+            time_string = BenchmarkTools.prettytime(time * 1e9)
+            println("$(neighborhood_searches_names[i])")
+            println("with $(join(sizes[iter], "x")) = $(prod(sizes[iter])) particles finished in $time_string\n")
+        end
+    end
+
+    plot(n_particles, times,
+         xaxis = :log, yaxis = :log,
+         xticks = (n_particles, n_particles),
+         xlabel = "#particles", ylabel = "Runtime [s]",
+         legend = :outerright, size = (750, 400), dpi = 200,
+         label = neighborhood_searches_names, title = title)
+end