Deactivate out of bounds particle (#1000)

* deactivate out of bounds particle

* check if update is necessary

* fix

* apply formatter

* fix tests

* fix tests again

* add unit tests

* implement suggestions

* fix tests

* add exact cell check

* add comments

* fix

* add comment

* change "this" to "the" in comment

* adapt comment

---------

Co-authored-by: LasNikas <niklas.nehe@web.de>

Author: LasNikas

src/general/neighborhood_search.jl

@@ -11,6 +11,44 @@ function PointNeighbors.foreach_point_neighbor(f, system, neighbor_system,
                            points, parallelization_backend)
 end
 
+deactivate_out_of_bounds_particles!(system, buffer, nhs, v, u, semi) = system
+
+function deactivate_out_of_bounds_particles!(system, buffer::SystemBuffer,
+                                             nhs::GridNeighborhoodSearch, v, u, semi)
+    deactivate_out_of_bounds_particles!(system, buffer, nhs, nhs.cell_list, v, u, semi)
+end
+
+function deactivate_out_of_bounds_particles!(system, buffer, nhs, cell_list, v, u, semi)
+    return system
+end
+
+# `GridNeighborhoodSearch` with a `FullGridCellList` requires a bounding box.
+# This function deactivates particles that move outside the bounding box to prevent
+# simulation crashes.
+# Note that deactivating particles is only possible in combination with a 'SystemBuffer'.
+function deactivate_out_of_bounds_particles!(system, ::SystemBuffer, nhs,
+                                             cell_list::FullGridCellList, v, u, semi)
+    @trixi_timeit timer() "deactivate out of bounds particle" begin
+        @threaded semi for particle in each_integrated_particle(system)
+            particle_position = current_coords(u, system, particle)
+            cell = PointNeighbors.cell_coords(particle_position, nhs)
+
+            # This is the same code as is used in PointNeighbors.jl in `check_cell_bounds`.
+            # It tests that particles are inside the inner grid (without the padding layer for neighbor query).
+            if !all(cell[i] in 2:(size(cell_list.linear_indices, i) - 1)
+                    for i in eachindex(cell))
+                deactivate_particle!(system, particle, v, u)
+            end
+        end
+
+        if count(system.buffer.active_particle) != system.buffer.active_particle_count[]
+            update_system_buffer!(system.buffer)
+        end
+    end
+
+    return system
+end
+
 @propagate_inbounds function foreach_neighbor(f, system_coords, neighbor_coords,
                                               neighborhood_search, backend, particle)
     PointNeighbors.foreach_neighbor(f, system_coords, neighbor_coords,

src/schemes/boundary/open_boundary/system.jl

@@ -312,6 +312,16 @@ end
     return du
 end
 
+function update_positions!(system::OpenBoundarySystem, v, u, v_ode, u_ode, semi, t)
+    nhs = get_neighborhood_search(system.fluid_system, system, semi)
+
+    # `GridNeighborhoodSearch` with a `FullGridCellList` requires a bounding box.
+    # This function deactivates particles that move outside the bounding box to prevent
+    # simulation crashes.
+    # Note that deactivating particles is only possible in combination with a 'SystemBuffer'.
+    deactivate_out_of_bounds_particles!(system, buffer(system), nhs, v, u, semi)
+end
+
 function update_boundary_interpolation!(system::OpenBoundarySystem, v, u, v_ode, u_ode,
                                         semi, t)
     update_boundary_model!(system, system.boundary_model, v, u, v_ode, u_ode, semi, t)

src/schemes/fluid/fluid.jl

@@ -115,6 +115,16 @@ end
 
 @inline acceleration_source(system::AbstractFluidSystem) = system.acceleration
 
+function update_positions!(system::AbstractFluidSystem, v, u, v_ode, u_ode, semi, t)
+    nhs = get_neighborhood_search(system, semi)
+
+    # `GridNeighborhoodSearch` with a `FullGridCellList` requires a bounding box.
+    # This function deactivates particles that move outside the bounding box to prevent
+    # simulation crashes.
+    # Note that deactivating particles is only possible in combination with a 'SystemBuffer'.
+    deactivate_out_of_bounds_particles!(system, buffer(system), nhs, v, u, semi)
+end
+
 function compute_density!(system, u, u_ode, semi, ::ContinuityDensity)
     # No density update with `ContinuityDensity`
     return system

test/general/general.jl

@@ -6,3 +6,4 @@ include("interpolation.jl")
 include("buffer.jl")
 include("util.jl")
 include("custom_quantities.jl")
+include("neighborhood_search.jl")

test/general/neighborhood_search.jl

@@ -0,0 +1,80 @@
+@testset verbose=true "Deactivate Out of Bounds Particles" begin
+    struct MockSystemOutOfBounds <: TrixiParticles.AbstractSystem{2}
+        buffer::TrixiParticles.SystemBuffer
+    end
+
+    TrixiParticles.nparticles(system::MockSystemOutOfBounds) = length(system.buffer.active_particle)
+    Base.eltype(system::MockSystemOutOfBounds) = Float64
+
+    @testset "Particles Inside Bounds" begin
+        # Setup: 5 particles, all inside bounds
+        buffer = TrixiParticles.SystemBuffer(5, 0)
+        system = MockSystemOutOfBounds(buffer)
+
+        u = [-0.5 0.0 0.5 -0.8 0.8
+             -0.5 -0.5 0.0 0.5 0.5]
+
+        cell_list = TrixiParticles.FullGridCellList(; min_corner=(-1.0, -1.0),
+                                                    max_corner=(1.0, 1.0),
+                                                    search_radius=0.1)
+        dummy_nhs = (; cell_size=0.1, periodic_box=nothing, cell_list)
+        semi = DummySemidiscretization()
+
+        # All particles should remain active
+        initial_count = count(buffer.active_particle)
+        TrixiParticles.deactivate_out_of_bounds_particles!(system, buffer, dummy_nhs,
+                                                           cell_list, u, u, semi)
+        @test count(buffer.active_particle) == initial_count
+    end
+
+    @testset "Particles Outside Bounds" begin
+        # Setup: 5 particles, some outside bounds
+        buffer = TrixiParticles.SystemBuffer(5, 0)
+        system = MockSystemOutOfBounds(buffer)
+
+        # Particles 3 and 5 are outside the bounds
+        u = [-0.5 0.0 2.0 -0.8 -2.0
+             -0.5 -0.5 0.0 0.5 0.5]
+
+        cell_list = TrixiParticles.FullGridCellList(; min_corner=(-1.0, -1.0),
+                                                    max_corner=(1.0, 1.0),
+                                                    search_radius=0.1)
+        dummy_nhs = (; cell_size=0.1, periodic_box=nothing, cell_list)
+        semi = DummySemidiscretization()
+
+        TrixiParticles.deactivate_out_of_bounds_particles!(system, buffer, dummy_nhs,
+                                                           cell_list, u, u, semi)
+
+        # Particles 3 and 5 should be deactivated
+        @test buffer.active_particle[3] == false
+        @test buffer.active_particle[5] == false
+        # Others should still be active
+        @test buffer.active_particle[1] == true
+        @test buffer.active_particle[2] == true
+        @test buffer.active_particle[4] == true
+
+        @test TrixiParticles.each_active_particle(system, buffer) == [1, 2, 4]
+    end
+
+    @testset "Edge Cases" begin
+        # Test the 1001//1000 padding logic
+        buffer = TrixiParticles.SystemBuffer(3, 0)
+        system = MockSystemOutOfBounds(buffer)
+
+        # Particles directly at the boundaries (should remain active)
+        u = [-1.0 1.0 0.0
+             -1.0 1.0 0.0]
+
+        cell_list = TrixiParticles.FullGridCellList(; min_corner=(-1.0, -1.0),
+                                                    max_corner=(1.0, 1.0),
+                                                    search_radius=0.1)
+        dummy_nhs = (; cell_size=0.1, periodic_box=nothing, cell_list)
+        semi = DummySemidiscretization()
+
+        TrixiParticles.deactivate_out_of_bounds_particles!(system, buffer, dummy_nhs,
+                                                           cell_list, u, u, semi)
+
+        # All should still be active
+        @test all(buffer.active_particle)
+    end
+end