diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl
index 8ac6d09bd..5fb40a762 100644
--- a/src/general/neighborhood_search.jl
+++ b/src/general/neighborhood_search.jl
@@ -15,8 +15,8 @@ function PointNeighbors.foreach_point_neighbor(f, system::GPUSystem, neighbor_sy
                                                neighborhood_search;
                                                points=each_moving_particle(system),
                                                parallel=true)
-    @threaded system for point in points
-        PointNeighbors.foreach_neighbor(f, system_coords, neighbor_coords,
-                                        neighborhood_search, point)
-    end
+    # For `GPUSystem`s, explicitly pass the backend, so a `GPUSystem` with a CPU
+    # backend will actually launch the KernelAbstractions.jl kernels on the CPU.
+    foreach_point_neighbor(f, system_coords, neighbor_coords, neighborhood_search;
+                           points, parallel=KernelAbstractions.get_backend(system_coords))
 end
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index ff17cc401..b1313dd44 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -644,20 +644,20 @@ function update_nhs!(neighborhood_search,
                      neighbor::Union{FluidSystem, TotalLagrangianSPHSystem},
                      u_system, u_neighbor)
     # The current coordinates of fluids and solids change over time
-    PointNeighbors.update!(neighborhood_search,
-                           current_coordinates(u_system, system),
-                           current_coordinates(u_neighbor, neighbor),
-                           points_moving=(true, true))
+    update!(neighborhood_search, system,
+            current_coordinates(u_system, system),
+            current_coordinates(u_neighbor, neighbor),
+            points_moving=(true, true))
 end
 
 function update_nhs!(neighborhood_search,
                      system::FluidSystem, neighbor::BoundarySPHSystem,
                      u_system, u_neighbor)
     # Boundary coordinates only change over time when `neighbor.ismoving[]`
-    PointNeighbors.update!(neighborhood_search,
-                           current_coordinates(u_system, system),
-                           current_coordinates(u_neighbor, neighbor),
-                           points_moving=(true, neighbor.ismoving[]))
+    update!(neighborhood_search, system,
+            current_coordinates(u_system, system),
+            current_coordinates(u_neighbor, neighbor),
+            points_moving=(true, neighbor.ismoving[]))
 end
 
 function update_nhs!(neighborhood_search,
@@ -667,10 +667,10 @@ function update_nhs!(neighborhood_search,
 
     # TODO: Update only `active_coordinates` of open boundaries.
     # Problem: Removing inactive particles from neighboring lists is necessary.
-    PointNeighbors.update!(neighborhood_search,
-                           current_coordinates(u_system, system),
-                           current_coordinates(u_neighbor, neighbor),
-                           points_moving=(true, true))
+    update!(neighborhood_search, system,
+            current_coordinates(u_system, system),
+            current_coordinates(u_neighbor, neighbor),
+            points_moving=(true, true))
 end
 
 function update_nhs!(neighborhood_search,
@@ -680,20 +680,20 @@ function update_nhs!(neighborhood_search,
 
     # TODO: Update only `active_coordinates` of open boundaries.
     # Problem: Removing inactive particles from neighboring lists is necessary.
-    PointNeighbors.update!(neighborhood_search,
-                           current_coordinates(u_system, system),
-                           current_coordinates(u_neighbor, neighbor),
-                           points_moving=(true, true))
+    update!(neighborhood_search, system,
+            current_coordinates(u_system, system),
+            current_coordinates(u_neighbor, neighbor),
+            points_moving=(true, true))
 end
 
 function update_nhs!(neighborhood_search,
                      system::TotalLagrangianSPHSystem, neighbor::FluidSystem,
                      u_system, u_neighbor)
     # The current coordinates of fluids and solids change over time
-    PointNeighbors.update!(neighborhood_search,
-                           current_coordinates(u_system, system),
-                           current_coordinates(u_neighbor, neighbor),
-                           points_moving=(true, true))
+    update!(neighborhood_search, system,
+            current_coordinates(u_system, system),
+            current_coordinates(u_neighbor, neighbor),
+            points_moving=(true, true))
 end
 
 function update_nhs!(neighborhood_search,
@@ -708,10 +708,10 @@ function update_nhs!(neighborhood_search,
                      u_system, u_neighbor)
     # The current coordinates of solids change over time.
     # Boundary coordinates only change over time when `neighbor.ismoving[]`.
-    PointNeighbors.update!(neighborhood_search,
-                           current_coordinates(u_system, system),
-                           current_coordinates(u_neighbor, neighbor),
-                           points_moving=(true, neighbor.ismoving[]))
+    update!(neighborhood_search, system,
+            current_coordinates(u_system, system),
+            current_coordinates(u_neighbor, neighbor),
+            points_moving=(true, neighbor.ismoving[]))
 end
 
 function update_nhs!(neighborhood_search,
@@ -734,10 +734,10 @@ function update_nhs!(neighborhood_search,
     #
     # Boundary coordinates only change over time when `neighbor.ismoving[]`.
     # The current coordinates of fluids and solids change over time.
-    PointNeighbors.update!(neighborhood_search,
-                           current_coordinates(u_system, system),
-                           current_coordinates(u_neighbor, neighbor),
-                           points_moving=(system.ismoving[], true))
+    update!(neighborhood_search, system,
+            current_coordinates(u_system, system),
+            current_coordinates(u_neighbor, neighbor),
+            points_moving=(system.ismoving[], true))
 end
 
 function update_nhs!(neighborhood_search,
@@ -746,30 +746,30 @@ function update_nhs!(neighborhood_search,
                      u_system, u_neighbor)
     # `system` coordinates only change over time when `system.ismoving[]`.
     # `neighbor` coordinates only change over time when `neighbor.ismoving[]`.
-    PointNeighbors.update!(neighborhood_search,
-                           current_coordinates(u_system, system),
-                           current_coordinates(u_neighbor, neighbor),
-                           points_moving=(system.ismoving[], neighbor.ismoving[]))
+    update!(neighborhood_search, system,
+            current_coordinates(u_system, system),
+            current_coordinates(u_neighbor, neighbor),
+            points_moving=(system.ismoving[], neighbor.ismoving[]))
 end
 
 function update_nhs!(neighborhood_search,
                      system::DEMSystem, neighbor::DEMSystem,
                      u_system, u_neighbor)
     # Both coordinates change over time
-    PointNeighbors.update!(neighborhood_search,
-                           current_coordinates(u_system, system),
-                           current_coordinates(u_neighbor, neighbor),
-                           points_moving=(true, true))
+    update!(neighborhood_search, system,
+            current_coordinates(u_system, system),
+            current_coordinates(u_neighbor, neighbor),
+            points_moving=(true, true))
 end
 
 function update_nhs!(neighborhood_search,
                      system::DEMSystem, neighbor::BoundaryDEMSystem,
                      u_system, u_neighbor)
     # DEM coordinates change over time, the boundary coordinates don't
-    PointNeighbors.update!(neighborhood_search,
-                           current_coordinates(u_system, system),
-                           current_coordinates(u_neighbor, neighbor),
-                           points_moving=(true, false))
+    update!(neighborhood_search, system,
+            current_coordinates(u_system, system),
+            current_coordinates(u_neighbor, neighbor),
+            points_moving=(true, false))
 end
 
 function update_nhs!(neighborhood_search,
@@ -788,6 +788,18 @@ function update_nhs!(neighborhood_search,
     return neighborhood_search
 end
 
+# Forward to PointNeighbors.jl
+function update!(neighborhood_search, system, x, y; points_moving=(true, false))
+    PointNeighbors.update!(neighborhood_search, x, y; points_moving)
+end
+
+# For `GPUSystem`s, explicitly pass the backend, so that a `GPUSystem` with a CPU
+# backend will actually launch the KernelAbstractions.jl kernels on the CPU.
+function update!(neighborhood_search, system::GPUSystem, x, y; points_moving=(true, false))
+    PointNeighbors.update!(neighborhood_search, x, y; points_moving,
+                           parallelization_backend=KernelAbstractions.get_backend(system))
+end
+
 function check_configuration(systems)
     foreach_system(systems) do system
         check_configuration(system, systems)