diff --git a/benchmarks/benchmarks.jl b/benchmarks/benchmarks.jl
index 854359f..6a16f4e 100644
--- a/benchmarks/benchmarks.jl
+++ b/benchmarks/benchmarks.jl
@@ -1,5 +1,6 @@
 include("count_neighbors.jl")
 include("n_body.jl")
+include("smoothed_particle_hydrodynamics.jl")
 include("update.jl")
 
 include("plot.jl")
diff --git a/benchmarks/smoothed_particle_hydrodynamics.jl b/benchmarks/smoothed_particle_hydrodynamics.jl
new file mode 100644
index 0000000..bb90c5a
--- /dev/null
+++ b/benchmarks/smoothed_particle_hydrodynamics.jl
@@ -0,0 +1,80 @@
+using PointNeighbors
+using TrixiParticles
+using BenchmarkTools
+
+"""
+    benchmark_wcsph(neighborhood_search, coordinates; parallel = true)
+
+A benchmark of the right-hand side of a full real-life Weakly Compressible
+Smoothed Particle Hydrodynamics (WCSPH) simulation with TrixiParticles.jl.
+This method is used to simulate an incompressible fluid.
+"""
+function benchmark_wcsph(neighborhood_search, coordinates; parallel = true)
+    density = 1000.0
+    fluid = InitialCondition(; coordinates, density, mass = 0.1)
+
+    # Compact support == smoothing length for the Wendland kernel
+    smoothing_length = PointNeighbors.search_radius(neighborhood_search)
+    if ndims(neighborhood_search) == 1
+        smoothing_kernel = SchoenbergCubicSplineKernel{1}()
+    else
+        smoothing_kernel = WendlandC2Kernel{ndims(neighborhood_search)}()
+    end
+
+    sound_speed = 10.0
+    state_equation = StateEquationCole(; sound_speed, reference_density = density,
+                                       exponent = 1)
+
+    fluid_density_calculator = ContinuityDensity()
+    viscosity = ArtificialViscosityMonaghan(alpha = 0.02, beta = 0.0)
+    density_diffusion = DensityDiffusionMolteniColagrossi(delta = 0.1)
+
+    fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator,
+                                               state_equation, smoothing_kernel,
+                                               smoothing_length, viscosity = viscosity,
+                                               density_diffusion = density_diffusion)
+
+    v = vcat(fluid.velocity, fluid.density')
+    u = copy(fluid.coordinates)
+    dv = zero(v)
+
+    # Initialize the system
+    TrixiParticles.initialize!(fluid_system, neighborhood_search)
+    TrixiParticles.compute_pressure!(fluid_system, v)
+
+    return @belapsed TrixiParticles.interact!($dv, $v, $u, $v, $u, $neighborhood_search,
+                                              $fluid_system, $fluid_system)
+end
+
+"""
+    benchmark_tlsph(neighborhood_search, coordinates; parallel = true)
+
+A benchmark of the right-hand side of a full real-life Total Lagrangian
+Smoothed Particle Hydrodynamics (TLSPH) simulation with TrixiParticles.jl.
+This method is used to simulate an elastic structure.
+"""
+function benchmark_tlsph(neighborhood_search, coordinates; parallel = true)
+    material = (density = 1000.0, E = 1.4e6, nu = 0.4)
+    solid = InitialCondition(; coordinates, density = material.density, mass = 0.1)
+
+    # Compact support == smoothing length for the Wendland kernel
+    smoothing_length = PointNeighbors.search_radius(neighborhood_search)
+    if ndims(neighborhood_search) == 1
+        smoothing_kernel = SchoenbergCubicSplineKernel{1}()
+    else
+        smoothing_kernel = WendlandC2Kernel{ndims(neighborhood_search)}()
+    end
+
+    solid_system = TotalLagrangianSPHSystem(solid, smoothing_kernel, smoothing_length,
+                                            material.E, material.nu)
+
+    v = copy(solid.velocity)
+    u = copy(solid.coordinates)
+    dv = zero(v)
+
+    # Initialize the system
+    TrixiParticles.initialize!(solid_system, neighborhood_search)
+
+    return @belapsed TrixiParticles.interact!($dv, $v, $u, $v, $u, $neighborhood_search,
+                                              $solid_system, $solid_system)
+end
diff --git a/test/Project.toml b/test/Project.toml
index f57fdf2..c826e5b 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -3,8 +3,10 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TrixiParticles = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 
 [compat]
 BenchmarkTools = "1"
 Plots = "1"
 Test = "1"
+TrixiParticles = "0.2"
diff --git a/test/benchmarks.jl b/test/benchmarks.jl
index 699e173..fe551b4 100644
--- a/test/benchmarks.jl
+++ b/test/benchmarks.jl
@@ -14,6 +14,14 @@
             @test_nowarn_mod plot_benchmarks(benchmark_n_body, size, 2)
         end
 
+        @testset verbose=true "`benchmark_wcsph`" begin
+            @test_nowarn_mod plot_benchmarks(benchmark_wcsph, size, 2)
+        end
+
+        @testset verbose=true "`benchmark_tlsph`" begin
+            @test_nowarn_mod plot_benchmarks(benchmark_tlsph, size, 2)
+        end
+
         @testset verbose=true "`benchmark_initialize`" begin
             @test_nowarn_mod plot_benchmarks(benchmark_initialize, size, 2)
         end