works on CPU but not on GPU

Author: LasNikas

src/general/buffer.jl

@@ -1,16 +1,24 @@
-struct SystemBuffer{AP, APC, EP}
-    active_particle       :: AP # Vector{Bool}
-    active_particle_count :: APC
-    eachparticle          :: EP # Vector{Int}
+struct SystemBuffer{VB, RI, VI}
+    active_particle       :: VB  # Vector{Bool}
+    active_particle_count :: RI  # Ref{Int}
+    candidates            :: VI  # Vector{Int}
+    particle_outside      :: VB  # Vector{Bool}
+    available_particles   :: VI  # Vector{Int}
+    next_particle         :: VI  # Vector{Int}
+    eachparticle          :: VI  # Vector{Int}
     buffer_size           :: Int
 end
 
 function SystemBuffer(active_size, buffer_size::Integer)
     # We cannot use a `BitVector` here, as writing to a `BitVector` is not thread-safe
     active_particle = vcat(fill(true, active_size), fill(false, buffer_size))
+    candidates = collect(eachindex(active_particle))
+    particle_outside = vcat(fill(false, active_size), fill(true, buffer_size))
+    available_particles = collect(eachindex(active_particle))
     eachparticle = collect(eachindex(active_particle))
 
-    return SystemBuffer(active_particle, Ref(active_size), eachparticle, buffer_size)
+    return SystemBuffer(active_particle, Ref(active_size), candidates, particle_outside,
+                        available_particles, Int[1], eachparticle, buffer_size)
 end
 
 allocate_buffer(initial_condition, ::Nothing) = initial_condition
@@ -42,7 +50,6 @@ end
 @inline function update_system_buffer!(buffer::SystemBuffer, semi)
     (; active_particle) = buffer
 
-
     buffer.active_particle_count[] = count(active_particle)
     buffer.eachparticle .= -1
 
@@ -62,31 +69,15 @@ end
     return buffer
 end
 
-@inline each_moving_particle(system, buffer) = view(buffer.eachparticle,
-                                                    1:buffer.active_particle_count[])
+@inline each_moving_particle(system,
+                             buffer) = view(buffer.eachparticle,
+                                            1:buffer.active_particle_count[])
 
 @inline active_coordinates(u, system, buffer) = view(u, :, buffer.active_particle)
 
-@inline active_particles(system, buffer) = view(buffer.eachparticle,
-                                                1:buffer.active_particle_count[])
-
-@inline function activate_next_particle(system)
-    (; active_particle) = system.buffer
-
-    for particle in eachindex(active_particle)
-        if !active_particle[particle]
-            # Activate this particle. The return value is the old value.
-            # If this is `true`, the particle was active before and we need to continue.
-            # This happens because a particle might have been activated by another thread
-            # between the condition and the line below.
-            was_active = PointNeighbors.Atomix.@atomicswap active_particle[particle] = true
-
-            !was_active && return particle
-        end
-    end
-
-    error("No buffer particles available")
-end
+@inline active_particles(system,
+                         buffer) = view(buffer.eachparticle,
+                                        1:buffer.active_particle_count[])
 
 @inline function deactivate_particle!(system, particle, u)
     (; active_particle) = system.buffer

src/schemes/boundary/open_boundary/method_of_characteristics.jl

@@ -25,7 +25,7 @@ end
 @inline function update_boundary_quantities!(system, boundary_model::BoundaryModelLastiwka,
                                              v, u, v_ode, u_ode, semi, t)
     (; density, pressure, cache, boundary_zone,
-    reference_velocity, reference_pressure, reference_density) = system
+     reference_velocity, reference_pressure, reference_density) = system
     (; flow_direction) = boundary_zone
 
     fluid_system = corresponding_fluid_system(system, semi)
@@ -155,7 +155,7 @@ end
 function evaluate_characteristics!(system, neighbor_system::FluidSystem,
                                    v, u, v_ode, u_ode, semi, t)
     (; volume, cache, boundary_zone, density, pressure,
-    reference_velocity, reference_pressure, reference_density) = system
+     reference_velocity, reference_pressure, reference_density) = system
     (; flow_direction) = boundary_zone
     (; characteristics) = cache
 

src/schemes/boundary/open_boundary/system.jl

@@ -281,31 +281,79 @@ function check_domain!(system, v, u, v_ode, u_ode, semi)
     u_fluid = wrap_u(u_ode, fluid_system, semi)
     v_fluid = wrap_v(v_ode, fluid_system, semi)
 
+    function map_sort!(a, b) # (a::Vector{Int}, b::Vector{Bool})
+        map!((true_value, idx) -> true_value ? -1 : idx, a, b, eachindex(b))
+
+        sort!(a, rev=true)
+
+        return a # [idx_1, idx_2, ..., -1, -1]
+    end
+
+    system.buffer.particle_outside .= false
+
     # Check the boundary particles whether they're leaving the boundary zone
     @threaded semi for particle in each_moving_particle(system)
         particle_coords = current_coords(u, system, particle)
 
         # Check if boundary particle is outside the boundary zone
         if !is_in_boundary_zone(boundary_zone, particle_coords)
-            convert_particle!(system, fluid_system, boundary_zone, particle,
-                              v, u, v_fluid, u_fluid)
+            system.buffer.particle_outside[particle] = true
         end
     end
 
+    # Basically: # candidates = [candidate_1, candidate_2, ..., -1, -1]
+    map_sort!(system.buffer.candidates, .!(system.buffer.particle_outside))
+
+    # Basically: `available_particles = [inactive_1, inactive_2, ..., -1, -1]`
+    map_sort!(fluid_system.buffer.available_particles, fluid_system.buffer.active_particle)
+
+    fluid_system.buffer.next_particle[1] = 0
+
+    # Copy buffer particle information to the fluid system
+    @threaded semi for buffer_id in 1:count(system.buffer.particle_outside)
+        particle = system.buffer.candidates[buffer_id]
+
+        fluid_id = PointNeighbors.Atomix.@atomic fluid_system.buffer.next_particle[1] += 1
+        new_particle = fluid_system.buffer.available_particles[fluid_id]
+
+        convert_particle!(system, fluid_system, boundary_zone, particle, new_particle,
+                          v, u, v_fluid, u_fluid)
+    end
+
     update_system_buffer!(system.buffer, semi)
     update_system_buffer!(fluid_system.buffer, semi)
 
+    fluid_system.buffer.particle_outside .= false
+
     # Check the fluid particles whether they're entering the boundary zone
     @threaded semi for fluid_particle in each_moving_particle(fluid_system)
         fluid_coords = current_coords(u_fluid, fluid_system, fluid_particle)
 
         # Check if fluid particle is in boundary zone
         if is_in_boundary_zone(boundary_zone, fluid_coords)
-            convert_particle!(fluid_system, system, boundary_zone, fluid_particle,
-                              v, u, v_fluid, u_fluid)
+            fluid_system.buffer.particle_outside[fluid_particle] = true
         end
     end
 
+    # Basically: # candidates = [candidate_1, candidate_2, ..., -1, -1]
+    map_sort!(fluid_system.buffer.candidates, .!(fluid_system.buffer.particle_outside))
+
+    # Basically: `available_particles = [inactive_1, inactive_2, ..., -1, -1]`
+    map_sort!(system.buffer.available_particles, system.buffer.active_particle)
+
+    system.buffer.next_particle[1] = 0
+
+    # Copy fluid particle information to the buffer
+    @threaded semi for fluid_id in 1:count(fluid_system.buffer.particle_outside)
+        particle = fluid_system.buffer.candidates[fluid_id]
+
+        buffer_id = PointNeighbors.Atomix.@atomic system.buffer.next_particle[1] += 1
+        new_particle = system.buffer.available_particles[buffer_id]
+
+        convert_particle!(fluid_system, system, boundary_zone, particle, new_particle,
+                          v, u, v_fluid, u_fluid)
+    end
+
     update_system_buffer!(system.buffer, semi)
     update_system_buffer!(fluid_system.buffer, semi)
 
@@ -314,21 +362,21 @@ end
 
 # Outflow particle is outside the boundary zone
 @inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system,
-                                   boundary_zone::BoundaryZone{OutFlow}, particle, v, u,
-                                   v_fluid, u_fluid)
+                                   boundary_zone::BoundaryZone{OutFlow},
+                                   particle, new_partcile, v, u, v_fluid, u_fluid)
     deactivate_particle!(system, particle, u)
 
     return system
 end
 
 # Inflow particle is outside the boundary zone
 @inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system,
-                                   boundary_zone::BoundaryZone{InFlow}, particle, v, u,
-                                   v_fluid, u_fluid)
+                                   boundary_zone::BoundaryZone{InFlow},
+                                   particle, new_partcile, v, u, v_fluid, u_fluid)
     (; spanning_set) = boundary_zone
 
     # Activate a new particle in simulation domain
-    transfer_particle!(fluid_system, system, particle, v_fluid, u_fluid, v, u)
+    transfer_particle!(fluid_system, system, particle, new_partcile, v_fluid, u_fluid, v, u)
 
     # Reset position of boundary particle
     for dim in 1:ndims(system)
@@ -340,8 +388,8 @@ end
 
 # Buffer particle is outside the boundary zone
 @inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system,
-                                   boundary_zone::BoundaryZone{BidirectionalFlow}, particle,
-                                   v, u, v_fluid, u_fluid)
+                                   boundary_zone::BoundaryZone{BidirectionalFlow},
+                                   particle, new_particle, v, u, v_fluid, u_fluid)
     relative_position = current_coords(u, system, particle) - boundary_zone.zone_origin
 
     # Check if particle is in- or outside the fluid domain.
@@ -353,7 +401,7 @@ end
     end
 
     # Activate a new particle in simulation domain
-    transfer_particle!(fluid_system, system, particle, v_fluid, u_fluid, v, u)
+    transfer_particle!(fluid_system, system, particle, new_particle, v_fluid, u_fluid, v, u)
 
     # Reset position of boundary particle
     for dim in 1:ndims(system)
@@ -365,19 +413,20 @@ end
 
 # Fluid particle is in boundary zone
 @inline function convert_particle!(fluid_system::FluidSystem, system,
-                                   boundary_zone, particle, v, u, v_fluid, u_fluid)
+                                   boundary_zone,
+                                   particle, new_particle, v, u, v_fluid, u_fluid)
     # Activate particle in boundary zone
-    transfer_particle!(system, fluid_system, particle, v, u, v_fluid, u_fluid)
+    transfer_particle!(system, fluid_system, particle, new_particle, v, u, v_fluid, u_fluid)
 
     # Deactivate particle in interior domain
     deactivate_particle!(fluid_system, particle, u_fluid)
 
     return fluid_system
 end
 
-@inline function transfer_particle!(system_new, system_old, particle_old,
+@inline function transfer_particle!(system_new, system_old, particle_old, particle_new,
                                     v_new, u_new, v_old, u_old)
-    particle_new = activate_next_particle(system_new)
+    system_new.buffer.active_particle[particle_new] = true
 
     # Transfer densities
     density = current_density(v_old, system_old, particle_old)