diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml
index 6a0a49bd..76e07ff0 100644
--- a/.github/workflows/Documenter.yml
+++ b/.github/workflows/Documenter.yml
@@ -53,6 +53,10 @@ jobs:
                     joinpath(output_directory, "machines_tlsph.png"))
             savefig(plot!(plot_machines_tlsph_deformation_grad(); title = ""),
                     joinpath(output_directory, "machines_tlsph_deformation_gradient.png"))
+            savefig(plot!(plot_implementations_wcsph(); title = ""),
+                    joinpath(output_directory, "implementations_wcsph.png"))
+            savefig(plot!(plot_update_strategies(); title = ""),
+                    joinpath(output_directory, "update_strategies.png"))
           '
       - name: Build and deploy
         env:
diff --git a/README.md b/README.md
index 3a8f291b..e968b4fa 100644
--- a/README.md
+++ b/README.md
@@ -21,8 +21,8 @@
 - Designed as a "playground" to easily switch between different implementations and data
   structures
 - Common API over all implementations
-- Extensive benchmark suite to study different implementations (work in progress)
-- GPU compatibility (work in progress)
+- Extensive benchmark suite to study different implementations
+- GPU compatibility
 
 | Implementation  | Description | Features | Query | Update | GPU-compatible |
 | ------------- | ------------- | --- | :--: | :--: | :--: |
@@ -32,16 +32,17 @@
 
 ## Benchmarks by Implementation
 
-The following benchmarks were conducted on an AMD Ryzen Threadripper 3990X using 128 threads.
+The following benchmarks were conducted on an Intel Xeon W9-3475X using 36 threads.
 
 Benchmark of a single force computation step of a Weakly Compressible SPH (WCSPH) simulation:
-![wcsph](https://github.com/trixi-framework/PointNeighbors.jl/assets/44124897/ad5c378b-9ce2-4e6f-91dc-1e0da379b91f)
-
-Benchmark of an incremental update similar to a WCSPH simulation (note the log scale):
-![update](https://github.com/trixi-framework/PointNeighbors.jl/assets/44124897/71eac5c9-6aa5-4267-bc0b-4057c89f8b12)
-
-Benchmark of a full right-hand side evaluation of a WCSPH simulation (note the log scale):
-![rhs](https://github.com/trixi-framework/PointNeighbors.jl/assets/44124897/ac328a96-1b9f-4319-a785-dce9d862fd70)
+![WCSPH benchmark by implementation](https://trixi-framework.org/PointNeighbors.jl/dev/assets/benchmarks/implementations_wcsph.png)
+
+Benchmark of different neighborhood search update strategies. The benchmark alternates
+between two perturbed point clouds. In 3D, only about 0.7% of the particles change cells
+between updates, which is similar to a dam break simulation.
+This is not relevant for `ParallelUpdate`, which reinitializes all particles
+in every update step.
+![Benchmark by update strategy](https://trixi-framework.org/PointNeighbors.jl/dev/assets/benchmarks/update_strategies.png)
 
 
 ## Benchmarks by Machine
diff --git a/benchmarks/plot_benchmarks.jl b/benchmarks/plot_benchmarks.jl
index 461e9159..02a41dfb 100644
--- a/benchmarks/plot_benchmarks.jl
+++ b/benchmarks/plot_benchmarks.jl
@@ -15,8 +15,9 @@ Keyword arguments are passed to `Plots.plot`. For example, use `title = "My titl
 
 # Examples
 ```julia
-include("benchmarks/benchmarks.jl")
-include("benchmarks/plot_benchmarks.jl")
+using PointNeighbors
+include(joinpath(pkgdir(PointNeighbors),  "benchmarks", "benchmarks.jl"));
+include(joinpath(pkgdir(PointNeighbors),  "benchmarks", "plot_benchmarks.jl"));
 
 n_particles_vec, times = run_benchmark_default(benchmark_count_neighbors, (10, 10), 3)
 plot_benchmark(n_particles_vec, times; title = "Count neighbors benchmark")
@@ -41,6 +42,8 @@ function plot_benchmark!(p, n_particles_vec, times; kwargs...)
 
     plot!(p, n_particles_vec, n_particles_vec ./ times .* 1e-6;
           xaxis = :log, xticks = (n_particles_vec, xticks), linewidth = 2,
+          # Make sure the plot starts at y = 0.
+          ylimits = (0, Inf), widen = true,
           xlabel = "#particles", ylabel = "million particles processed per second",
           legend = :outerright, size = (700, 350), dpi = 600, margin = 4 * Plots.mm,
           palette = palette(:tab10), kwargs...)
@@ -50,6 +53,7 @@ end
 # Run these benchmarks like this. The type of `search_radius_factor` determines
 # if the benchmark is run in Float32 or Float64.
 #
+# using PointNeighbors
 # include(joinpath(pkgdir(PointNeighbors),  "benchmarks", "benchmarks.jl"));
 # _, times = run_benchmark_full_grid(benchmark_wcsph, (10, 10, 10), 9,
 #                                    search_radius_factor=3.0f0,
@@ -86,7 +90,22 @@ benchmark_runtimes = (n_particles = [
                       n_particles_tlsph_deformation_grad_5090 = [
                           1000, 4096, 15625, 64000, 250047, 1030301, 4096000, 16387064,
                           42875000],
-                      tlsph_deformation_grad_5090_fp32 = [5.6569e-5; 6.894e-5; 7.6672e-5; 8.3602e-5; 0.000313116; 0.00099187; 0.003671176; 0.015459844; 0.041021319;;])
+                      tlsph_deformation_grad_5090_fp32 = [5.6569e-5; 6.894e-5; 7.6672e-5; 8.3602e-5; 0.000313116; 0.00099187; 0.003671176; 0.015459844; 0.041021319;;],
+                      # Benchmarking the difference between implementations on an
+                      # Intel Xeon W9-3475X (x36).
+                      # _, times = run_benchmark_default(benchmark_wcsph, ...)
+                      wcsph_w9_3475x_dictionary = [9.539e-5; 0.000507708; 0.001999838; 0.008442291; 0.034369512; 0.145935124; 0.581081578; 2.341472385; 9.445755411;;],
+                      wcsph_w9_3475x_fullgrid = [6.0765e-5; 0.000402875; 0.001569317; 0.006578832; 0.026040483; 0.108027132; 0.434874618; 1.755279267; 7.07428452;;],
+                      wcsph_w9_3475x_precomputed = [3.7762e-5; 0.000148563; 0.00054964; 0.002274807; 0.00899912; 0.037950566; 0.151298558; 0.607103175; 2.440462785;;],
+                      # NaN values are just placeholders because the benchmark took too long.
+                      wcsph_w9_3475x_trivial = [6.3518e-5; 0.000669186; 0.008641863; 0.135930064; 2.468379663; 43.607518925; NaN; NaN; NaN;;],
+                      # Benchmarking the difference between update strategies on an
+                      # Intel Xeon W9-3475X (x36).
+                      # _, times = run_benchmark_updates((10, 10, 10), 9)
+                      update_w9_3475x_parallel = [1.42395e-5; 1.9205e-5; 3.3569e-5; 7.6216e-5; 0.0001972665; 0.000657154; 0.002366056; 0.008558718; 0.0328478125;;],
+                      update_w9_3475x_parallel_incremental = [1.05465e-5; 1.41955e-5; 2.36415e-5; 6.3423e-5; 0.000185218; 0.0007833225; 0.0052182165; 0.025595217; 0.103340834;;],
+                      update_w9_3475x_semi_parallel = [1.00805e-5; 1.78505e-5; 4.28215e-5; 0.000153458; 0.0006931965; 0.0026806445; 0.0158252845; 0.079710013; 0.329653309;;],
+                      update_w9_3475x_precomputed = [0.000156829; 0.0005617385; 0.0021646195; 0.009024255; 0.0354495175; 0.1443668715; 0.5821228515; 2.341333004; 9.381562997;;])
 
 function plot_machines_wcsph()
     times = hcat(benchmark_runtimes.wcsph_5090_fp32,
@@ -132,3 +151,33 @@ function plot_machines_tlsph_deformation_grad()
 
     plot_benchmark!(p, benchmark_runtimes.n_particles, times; label = names)
 end
+
+function plot_implementations_wcsph()
+    times = hcat(benchmark_runtimes.wcsph_w9_3475x_precomputed,
+                 benchmark_runtimes.wcsph_w9_3475x_fullgrid,
+                 benchmark_runtimes.wcsph_w9_3475x_dictionary,
+                 benchmark_runtimes.wcsph_w9_3475x_trivial)
+
+    names = ["PrecomputedNeighborhoodSearch";;
+             "GNHS & FullGridCellList";;
+             "GridNeighborhoodSearch";;
+             "TrivialNeighborhoodSearch";;]
+
+    plot_benchmark(benchmark_runtimes.n_particles, times; label = names,
+                   title = "WCSPH on Intel Xeon W9-3475X (x36)")
+end
+
+function plot_update_strategies()
+    times = hcat(benchmark_runtimes.update_w9_3475x_parallel,
+                 benchmark_runtimes.update_w9_3475x_parallel_incremental,
+                 benchmark_runtimes.update_w9_3475x_semi_parallel,
+                 benchmark_runtimes.update_w9_3475x_precomputed)
+
+    names = ["GNHS & ParallelUpdate";;
+             "GNHS & ParallelIncrementalUpdate";;
+             "GNHS & SemiParallelUpdate";;
+             "PrecomputedNeighborhoodSearch";;]
+
+    plot_benchmark(benchmark_runtimes.n_particles, times; label = names,
+                   title = "Update strategies on Intel Xeon W9-3475X (x36)")
+end
diff --git a/benchmarks/run_benchmarks.jl b/benchmarks/run_benchmarks.jl
index 8e671e0d..11f5ee79 100644
--- a/benchmarks/run_benchmarks.jl
+++ b/benchmarks/run_benchmarks.jl
@@ -38,7 +38,8 @@ See also
                             of `search_radius_factor` determines if the benchmark
                             is run in single or double precision.
 - `parallelization_backend = PolyesterBackend()`: Parallelization strategy to use. See
-                            [`@threaded`](@ref) for a list of available backends.
+                            [`PointNeighbors.@threaded`](@ref) for a list of available
+                            backends.
 - `names = ["Neighborhood search 1" ...]`: Names of the neighborhood searches used in the
                             benchmark output.
 - `seed = 1`:               Seed to perturb the point positions. Different seeds yield
@@ -51,7 +52,8 @@ See also
 
 # Examples
 ```julia
-include("benchmarks/benchmarks.jl")
+using PointNeighbors
+include(joinpath(pkgdir(PointNeighbors),  "benchmarks", "benchmarks.jl"));
 
 run_benchmark(benchmark_count_neighbors, (10, 10), 3,
               [TrivialNeighborhoodSearch{2}(), GridNeighborhoodSearch{2}()])
@@ -113,7 +115,8 @@ function run_benchmark(benchmark, n_points_per_dimension, iterations, neighborho
 end
 
 """
-    run_benchmark_default(benchmark, n_points_per_dimension, iterations; kwargs...)
+    run_benchmark_default(benchmark, n_points_per_dimension, iterations;
+                          max_neighbors = 128, kwargs...)
 
 Shortcut to call [`run_benchmark`](@ref) with the most commonly used neighborhood search
 implementations:
@@ -132,16 +135,23 @@ implementations:
 - `iterations`:             Number of refinement iterations
 
 # Keywords
-See [`run_benchmark`](@ref) for a list of available keywords.
+- `max_neighbors = 128`: Maximum neighbor-list capacity for the precomputed neighborhood
+                         search. This needs to be increased for a larger
+                         `search_radius_factor` and can be decreased if the benchmark runs
+                         out of memory.
+
+See [`run_benchmark`](@ref) for a list of additional available keywords.
 
 # Examples
 ```julia
-include("benchmarks/benchmarks.jl")
+using PointNeighbors
+include(joinpath(pkgdir(PointNeighbors),  "benchmarks", "benchmarks.jl"));
 
 run_benchmark_default(benchmark_n_body, (10, 10), 3)
 ```
 """
-function run_benchmark_default(benchmark, n_points_per_dimension, iterations; kwargs...)
+function run_benchmark_default(benchmark, n_points_per_dimension, iterations;
+                               max_neighbors = 128, kwargs...)
     NDIMS = length(n_points_per_dimension)
     min_corner = 0.0f0 .* n_points_per_dimension
     max_corner = Float32.(n_points_per_dimension ./ maximum(n_points_per_dimension))
@@ -151,7 +161,7 @@ function run_benchmark_default(benchmark, n_points_per_dimension, iterations; kw
         GridNeighborhoodSearch{NDIMS}(search_radius = 0.0f0,
                                       cell_list = FullGridCellList(; search_radius = 0.0f0,
                                                                    min_corner, max_corner)),
-        PrecomputedNeighborhoodSearch{NDIMS}()
+        PrecomputedNeighborhoodSearch{NDIMS}(; max_neighbors, sort_neighbor_lists = true)
     ]
 
     names = ["GridNeighborhoodSearch";;
@@ -163,7 +173,8 @@ function run_benchmark_default(benchmark, n_points_per_dimension, iterations; kw
 end
 
 """
-    run_benchmark_gpu(benchmark, n_points_per_dimension, iterations; kwargs...)
+    run_benchmark_gpu(benchmark, n_points_per_dimension, iterations;
+                      max_neighbors = 128, kwargs...)
 
 Shortcut to call [`run_benchmark`](@ref) with all GPU-compatible neighborhood search
 implementations:
@@ -181,20 +192,28 @@ implementations:
 - `iterations`:             Number of refinement iterations
 
 # Keywords
-See [`run_benchmark`](@ref) for a list of available keywords.
+- `max_neighbors = 128`: Maximum neighbor-list capacity for the precomputed neighborhood
+                         search. This needs to be increased for a larger
+                         `search_radius_factor` and can be decreased if the benchmark runs
+                         out of memory.
+
+See [`run_benchmark`](@ref) for a list of additional available keywords.
 
 # Examples
 ```julia
-include("benchmarks/benchmarks.jl")
+using PointNeighbors
+include(joinpath(pkgdir(PointNeighbors),  "benchmarks", "benchmarks.jl"));
 
 run_benchmark_gpu(benchmark_n_body, (10, 10), 3)
 ```
 """
 function run_benchmark_gpu(benchmark, n_points_per_dimension, iterations;
-                           parallelization_backend = PolyesterBackend(), kwargs...)
+                           parallelization_backend = PolyesterBackend(),
+                           max_neighbors = 128, kwargs...)
     grid_nhs = create_full_grid_neighborhood_search(n_points_per_dimension)
     precomputed_nhs = create_precomputed_neighborhood_search(grid_nhs,
-                                                             parallelization_backend)
+                                                             parallelization_backend;
+                                                             max_neighbors)
     neighborhood_searches = (grid_nhs, precomputed_nhs)
 
     names = ["GridNeighborhoodSearch with FullGridCellList";;
@@ -227,7 +246,8 @@ See [`run_benchmark`](@ref) for a list of available keywords.
 
 # Examples
 ```julia
-include("benchmarks/benchmarks.jl")
+using PointNeighbors
+include(joinpath(pkgdir(PointNeighbors),  "benchmarks", "benchmarks.jl"));
 
 run_benchmark_full_grid(benchmark_n_body, (10, 10), 3)
 ```
@@ -243,7 +263,8 @@ function run_benchmark_full_grid(benchmark, n_points_per_dimension, iterations;
 end
 
 """
-    run_benchmark_precomputed(benchmark, n_points_per_dimension, iterations; kwargs...)
+    run_benchmark_precomputed(benchmark, n_points_per_dimension, iterations;
+                              max_neighbors = 128, kwargs...)
 
 Shortcut to call [`run_benchmark`](@ref) with a `PrecomputedNeighborhoodSearch`.
 This is the neighborhood search implementation that is used in TrixiParticles.jl for
@@ -261,20 +282,28 @@ Use this function to benchmark and profile TrixiParticles.jl kernels.
 - `iterations`:             Number of refinement iterations
 
 # Keywords
-See [`run_benchmark`](@ref) for a list of available keywords.
+- `max_neighbors = 128`: Maximum neighbor-list capacity for the precomputed neighborhood
+                         search. This needs to be increased for a larger
+                         `search_radius_factor` and can be decreased if the benchmark runs
+                         out of memory.
+
+See [`run_benchmark`](@ref) for a list of additional available keywords.
 
 # Examples
 ```julia
-include("benchmarks/benchmarks.jl")
+using PointNeighbors
+include(joinpath(pkgdir(PointNeighbors),  "benchmarks", "benchmarks.jl"));
 
 run_benchmark_precomputed(benchmark_n_body, (10, 10), 3)
 ```
 """
 function run_benchmark_precomputed(benchmark, n_points_per_dimension, iterations;
-                                   parallelization_backend = PolyesterBackend(), kwargs...)
+                                   parallelization_backend = PolyesterBackend(),
+                                   max_neighbors = 128, kwargs...)
     grid_nhs = create_full_grid_neighborhood_search(n_points_per_dimension)
     precomputed_nhs = create_precomputed_neighborhood_search(grid_nhs,
-                                                             parallelization_backend)
+                                                             parallelization_backend;
+                                                             max_neighbors)
     neighborhood_searches = (precomputed_nhs,)
 
     names = ["PrecomputedNeighborhoodSearch";;]
@@ -283,21 +312,85 @@ function run_benchmark_precomputed(benchmark, n_points_per_dimension, iterations
                   neighborhood_searches; names, parallelization_backend, kwargs...)
 end
 
-function create_full_grid_neighborhood_search(n_points_per_dimension)
+"""
+    run_benchmark_updates(n_points_per_dimension, iterations;
+                          max_neighbors = 128, kwargs...)
+
+Benchmark [`benchmark_update_alternating`](@ref) with the update strategies
+[`ParallelUpdate`](@ref), [`ParallelIncrementalUpdate`](@ref), and
+[`SemiParallelUpdate`](@ref), as well as with a [`PrecomputedNeighborhoodSearch`](@ref)
+using [`ParallelUpdate`](@ref) internally.
+
+The benchmark alternates between two perturbed point clouds. In 3D, approximately 0.7% of
+the particles change cells between them, representative of an SPH update. This fraction is
+not relevant for [`ParallelUpdate`](@ref), which reinitializes all particles on every
+update.
+
+Returns `(n_particles_vec, times)` as described for [`run_benchmark`](@ref).
+
+# Arguments
+- `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
+- `max_neighbors = 128`: Maximum neighbor-list capacity for the precomputed neighborhood
+                         search. This needs to be increased for a larger
+                         `search_radius_factor` and can be decreased if the benchmark runs
+                         out of memory.
+
+See [`run_benchmark`](@ref) for a list of additional available keywords.
+
+# Examples
+```julia
+using PointNeighbors
+include(joinpath(pkgdir(PointNeighbors),  "benchmarks", "benchmarks.jl"));
+
+run_benchmark_updates((10, 10, 10), 3)
+```
+"""
+function run_benchmark_updates(n_points_per_dimension, iterations;
+                               parallelization_backend = PolyesterBackend(),
+                               max_neighbors = 128, kwargs...)
+    parallel = create_full_grid_neighborhood_search(n_points_per_dimension;
+                                                    update_strategy = ParallelUpdate())
+    parallel_incremental = create_full_grid_neighborhood_search(n_points_per_dimension;
+                                                                update_strategy = ParallelIncrementalUpdate())
+    semi_parallel = create_full_grid_neighborhood_search(n_points_per_dimension;
+                                                         update_strategy = SemiParallelUpdate())
+    precomputed_nhs = create_precomputed_neighborhood_search(parallel,
+                                                             parallelization_backend;
+                                                             max_neighbors)
+    neighborhood_searches = (parallel, parallel_incremental, semi_parallel, precomputed_nhs)
+
+    names = ["GNHS with ParallelUpdate";;
+             "GNHS with ParallelIncrementalUpdate";;
+             "GNHS with SemiParallelUpdate";;
+             "PrecomputedNeighborhoodSearch";;]
+
+    run_benchmark(benchmark_update_alternating, n_points_per_dimension, iterations,
+                  neighborhood_searches; names, parallelization_backend, kwargs...)
+end
+
+function create_full_grid_neighborhood_search(n_points_per_dimension;
+                                              update_strategy = ParallelUpdate())
     NDIMS = length(n_points_per_dimension)
 
     min_corner = 0.0f0 .* n_points_per_dimension
     max_corner = Float32.(n_points_per_dimension ./ maximum(n_points_per_dimension))
     cell_list = FullGridCellList(; search_radius = 0.0f0, min_corner, max_corner)
     return GridNeighborhoodSearch{NDIMS}(; search_radius = 0.0f0, cell_list,
-                                         update_strategy = ParallelUpdate())
+                                         update_strategy)
 end
 
-function create_precomputed_neighborhood_search(grid_nhs, parallelization_backend)
+function create_precomputed_neighborhood_search(grid_nhs, parallelization_backend;
+                                                max_neighbors = 128)
     NDIMS = ndims(grid_nhs)
     transpose_backend = parallelization_backend isa PointNeighbors.KernelAbstractions.GPU
     return PrecomputedNeighborhoodSearch{NDIMS}(; search_radius = 0.0f0,
-                                                max_neighbors = 128,
+                                                max_neighbors,
                                                 update_neighborhood_search = grid_nhs,
+                                                sort_neighbor_lists = true,
                                                 transpose_backend)
 end
diff --git a/benchmarks/smoothed_particle_hydrodynamics.jl b/benchmarks/smoothed_particle_hydrodynamics.jl
index 8019c3c0..6275bbb3 100644
--- a/benchmarks/smoothed_particle_hydrodynamics.jl
+++ b/benchmarks/smoothed_particle_hydrodynamics.jl
@@ -55,9 +55,13 @@ function benchmark_wcsph(neighborhood_search, coordinates;
                              mass = convert(ELTYPE, 0.1) * particle_spacing,
                              particle_spacing)
 
+    # Perturb the initial velocity so that approaching particle pairs exercise the
+    # `vr < 0` branch of ArtificialViscosityMonaghan.
+    fluid.velocity .+= convert(ELTYPE, 1.0e-3) .* randn(ELTYPE, size(fluid.velocity))
+
     # Make sure that the computed forces are not all zero
     for i in eachindex(fluid.density)
-        fluid.density[i] += rand(eltype(fluid.density))
+        fluid.density[i] += randn(eltype(fluid.density))
     end
 
     sound_speed = convert(ELTYPE, 10.0)
diff --git a/benchmarks/update.jl b/benchmarks/update.jl
index 79fb8bb2..8fec789b 100644
--- a/benchmarks/update.jl
+++ b/benchmarks/update.jl
@@ -24,16 +24,21 @@ A very simple benchmark for neighborhood search update, alternating between two
 perturbed point clouds.
 
 This is a good benchmark for incremental updates, since most particles stay in their cells.
+Similar to 3D dam break simulations, ~0.5% of the particles change their cell during
+an update in 2D and ~0.7% in 3D.
 """
 function benchmark_update_alternating(neighborhood_search, coordinates;
                                       parallelization_backend = default_backend(coordinates))
     coordinates2 = copy(coordinates)
-    # Perturb all coordinates with a perturbation factor of `0.015`.
+    # Perturb all coordinates with a perturbation factor of `4e-4`.
     # This factor was tuned so that ~0.5% of the particles change their cell during an
     # update in 2D and ~0.7% in 3D.
     # These values are the same as the experimentally computed averages in 2D and 3D SPH
     # dam break simulations. So this benchmark replicates a real-life SPH update.
-    perturb!(coordinates2, 4e-4 * PointNeighbors.search_radius(neighborhood_search))
+    # Use the same perturbation every time to make comparisons fair and reproducible.
+    rng = Random.Xoshiro(1)
+    perturb!(coordinates2,
+             4e-4 * PointNeighbors.search_radius(neighborhood_search); rng)
 
     function update_alternating!(neighborhood_search, coordinates, coordinates2,
                                  parallelization_backend)
diff --git a/docs/Project.toml b/docs/Project.toml
index 2be4453f..02cf2781 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,13 +1,17 @@
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+TrixiParticles = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 
 [compat]
 Adapt = "4"
+BenchmarkTools = "1"
 Documenter = "1"
 KernelAbstractions = "0.9"
 Literate = "2"
 Plots = "1"
+TrixiParticles = "0.5"
diff --git a/docs/make.jl b/docs/make.jl
index 4a24604f..b25bc322 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -12,6 +12,11 @@ end
 
 using PointNeighbors
 
+# The benchmark files are scripts that users include into `Main`. Load them the same way
+# here so their docstrings are documented under their actual `Main.<name>` bindings.
+include(joinpath(trixibase_root_dir, "benchmarks", "benchmarks.jl"))
+include(joinpath(trixibase_root_dir, "benchmarks", "plot_benchmarks.jl"))
+
 # Define module-wide setups such that the respective modules are available in doctests
 DocMeta.setdocmeta!(PointNeighbors, :DocTestSetup,
                     :(using PointNeighbors); recursive = true)
@@ -76,7 +81,8 @@ makedocs(modules = [PointNeighbors],
                  "GPU Usage" => joinpath("tutorials", "tut_gpu_usage.md"),
                  "Advanced Usage" => joinpath("tutorials", "tut_advanced_usage.md")
              ],
-             "API reference" => "reference.md",
+             "API Reference" => "reference.md",
+             "Benchmark Suite" => "benchmarks.md",
              "Authors" => "authors.md",
              "License" => "license.md"
          ])
diff --git a/docs/src/benchmarks.md b/docs/src/benchmarks.md
new file mode 100644
index 00000000..374b564a
--- /dev/null
+++ b/docs/src/benchmarks.md
@@ -0,0 +1,45 @@
+# Benchmark Suite
+
+```@meta
+CurrentModule = Main
+```
+
+PointNeighbors.jl includes a benchmark suite for comparing neighborhood search
+implementations and update strategies. Load it with
+
+```julia
+using PointNeighbors, BenchmarkTools, TrixiParticles
+include(joinpath(pkgdir(PointNeighbors),  "benchmarks", "benchmarks.jl"));
+```
+
+## Benchmark Runners
+
+```@autodocs
+Modules = [Main]
+Pages = ["run_benchmarks.jl"]
+Order = [:function]
+```
+
+## Benchmark Workloads
+
+```@autodocs
+Modules = [Main]
+Pages = ["count_neighbors.jl", "n_body.jl",
+         "smoothed_particle_hydrodynamics.jl", "update.jl"]
+Order = [:function]
+```
+
+## Plotting
+
+Load the plotting utilities with
+
+```julia
+using PointNeighbors, Plots
+include(joinpath(pkgdir(PointNeighbors), "benchmarks", "plot_benchmarks.jl"));
+```
+
+```@autodocs
+Modules = [Main]
+Pages = ["plot_benchmarks.jl"]
+Order = [:function]
+```
diff --git a/test/benchmarks.jl b/test/benchmarks.jl
index ccae6fc5..4afeb6ef 100644
--- a/test/benchmarks.jl
+++ b/test/benchmarks.jl
@@ -7,10 +7,13 @@
 
     @testset verbose=true "$(length(size))D" for size in [(50,), (10, 10), (5, 5, 5)]
         @testset verbose=true "`benchmark_count_neighbors`" begin
-            @trixi_test_nowarn run_benchmark_default(benchmark_count_neighbors, size, 2)
-            @trixi_test_nowarn run_benchmark_gpu(benchmark_count_neighbors, size, 2)
+            @trixi_test_nowarn run_benchmark_default(benchmark_count_neighbors, size, 2;
+                                                     max_neighbors = 128)
+            @trixi_test_nowarn run_benchmark_gpu(benchmark_count_neighbors, size, 2;
+                                                 max_neighbors = 128)
             @trixi_test_nowarn run_benchmark_full_grid(benchmark_count_neighbors, size, 2)
-            @trixi_test_nowarn run_benchmark_precomputed(benchmark_count_neighbors, size, 2)
+            @trixi_test_nowarn run_benchmark_precomputed(benchmark_count_neighbors, size,
+                                                         2; max_neighbors = 128)
         end
 
         @testset verbose=true "`benchmark_n_body`" begin
@@ -51,4 +54,8 @@
             @trixi_test_nowarn run_benchmark_gpu(benchmark_update_alternating, size, 2)
         end
     end
+
+    @testset verbose=true "`run_benchmark_updates`" begin
+        @trixi_test_nowarn run_benchmark_updates((5, 5, 5), 2; max_neighbors = 128)
+    end
 end;
diff --git a/test/point_cloud.jl b/test/point_cloud.jl
index 4f8a44f0..b0af830b 100644
--- a/test/point_cloud.jl
+++ b/test/point_cloud.jl
@@ -5,7 +5,7 @@ function point_cloud(n_points_per_dimension, search_radius;
                      seed = 1, perturbation_factor_position = 1.0,
                      shuffle = false, sort = !shuffle)
     # Fixed seed to ensure reproducibility
-    Random.seed!(seed)
+    rng = Random.Xoshiro(seed)
 
     n_dims = length(n_points_per_dimension)
     coordinates = Array{Float64}(undef, n_dims, prod(n_points_per_dimension))
@@ -23,20 +23,14 @@ function point_cloud(n_points_per_dimension, search_radius;
         # grid, 1.6 for a 600 x 200 grid and 1.26 for a 1200 x 400 grid.
         # This is consistent with the standard deviation in a vortex street simulation.
         # The benchmark results are also consistent with the timer output of the simulation.
-        point_coords += perturbation_factor_position * 0.05 * randn(typeof(point_coords))
+        point_coords += perturbation_factor_position * 0.05 *
+                        randn(rng, typeof(point_coords))
 
         coordinates[:, i] .= point_coords
         cell_coords[i] = PointNeighbors.nonperiodic_cell_coords(point_coords, nothing,
                                                                 cell_size) .+ 1
     end
 
-    # A standard deviation of 0.05 in the particle coordinates
-    # corresponds to a standard deviation of 2 in the number of neighbors for a 300 x 100
-    # grid, 1.6 for a 600 x 200 grid and 1.26 for a 1200 x 400 grid.
-    # This is consistent with the standard deviation in a vortex street simulation.
-    # The benchmark results are also consistent with the timer output of the simulation.
-    perturb!(coordinates, perturbation_factor_position * 0.05)
-
     # Sort by the cell coordinates of the perturbed points.
     if sort
         if shuffle
@@ -50,16 +44,16 @@ function point_cloud(n_points_per_dimension, search_radius;
         coordinates .= coordinates[:, permutation]
     elseif shuffle
         # Sort randomly
-        permutation = Random.shuffle(axes(coordinates, 2))
+        permutation = Random.shuffle(rng, axes(coordinates, 2))
         coordinates .= coordinates[:, permutation]
     end
 
     return coordinates
 end
 
-function perturb!(data, std_deviation)
+function perturb!(data, std_deviation; rng = Random.default_rng())
     for i in eachindex(data)
-        data[i] += std_deviation * randn()
+        data[i] += std_deviation * randn(rng)
     end
 
     return data