Implement set operations for initial conditions (#202)

* Add serial `for_particle_neighbor` * Implement set operations * Fix dummy particles test * Add tests * Add vortex street example * Fix example * Remove vortex street example --------- Co-authored-by: LasNikas <[email protected]>
trixi-framework · Aug 16, 2023 · aaeda28 · aaeda28
1 parent d998f8a
commit aaeda28
Show file tree

Hide file tree

Showing 11 changed files with 378 additions and 79 deletions.
diff --git a/examples/fsi/dam_break_2d.jl b/examples/fsi/dam_break_2d.jl
@@ -80,7 +80,7 @@ fixed_particles = RectangularShape(solid_particle_spacing,
                                    (2water_width, 0.0),
                                    solid_density, tlsph=true)
 
-solid = InitialCondition(plate, fixed_particles)
+solid = union(plate, fixed_particles)
 
 # ==========================================================================================
 # ==== Boundary models

diff --git a/examples/fsi/dam_break_gate_2d.jl b/examples/fsi/dam_break_gate_2d.jl
@@ -93,7 +93,7 @@ fixed_particles = RectangularShape(solid_particle_spacing,
                                    (n_particles_x, 1), (plate_position, 0.0),
                                    solid_density, tlsph=true)
 
-solid = InitialCondition(plate, fixed_particles)
+solid = union(plate, fixed_particles)
 
 # ==========================================================================================
 # ==== Boundary models

diff --git a/examples/fsi/falling_water_column_2d.jl b/examples/fsi/falling_water_column_2d.jl
@@ -65,7 +65,7 @@ n_particles_per_dimension = (round(Int, length_beam / solid_particle_spacing) +
 beam = RectangularShape(solid_particle_spacing, n_particles_per_dimension,
                         (0.0, 0.0), solid_density, tlsph=true)
 
-solid = InitialCondition(beam, fixed_particles)
+solid = union(beam, fixed_particles)
 
 # ==========================================================================================
 # ==== Boundary models

diff --git a/examples/solid/oscillating_beam_2d.jl b/examples/solid/oscillating_beam_2d.jl
@@ -41,7 +41,7 @@ n_particles_per_dimension = (round(Int, length_beam / particle_spacing) +
 beam = RectangularShape(particle_spacing, n_particles_per_dimension,
                         (0.0, 0.0), particle_density, tlsph=true)
 
-solid = InitialCondition(beam, fixed_particles)
+solid = union(beam, fixed_particles)
 
 # ==========================================================================================
 # ==== Systems

diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl
@@ -1,10 +1,12 @@
 struct InitialCondition{ELTYPE}
-    coordinates :: Array{ELTYPE, 2}
-    velocity    :: Array{ELTYPE, 2}
-    mass        :: Array{ELTYPE, 1}
-    density     :: Array{ELTYPE, 1}
+    particle_spacing :: ELTYPE
+    coordinates      :: Array{ELTYPE, 2}
+    velocity         :: Array{ELTYPE, 2}
+    mass             :: Array{ELTYPE, 1}
+    density          :: Array{ELTYPE, 1}
 
-    function InitialCondition(coordinates, velocities, masses, densities)
+    function InitialCondition(coordinates, velocities, masses, densities;
+                              particle_spacing=-1.0)
         if size(coordinates) != size(velocities)
             throw(ArgumentError("`coordinates` and `velocities` must be of the same size"))
         end
@@ -14,22 +16,8 @@ struct InitialCondition{ELTYPE}
                                 "`size(coordinates, 2) == length(masses) == length(densities)`"))
         end
 
-        return new{eltype(coordinates)}(coordinates, velocities, masses, densities)
-    end
-
-    function InitialCondition(initial_conditions...)
-        NDIMS = size(first(initial_conditions).coordinates, 1)
-        if any(ic -> size(ic.coordinates, 1) != NDIMS, initial_conditions)
-            throw(ArgumentError("all passed initial conditions must have the same dimensionality"))
-        end
-
-        coordinates = hcat((ic.coordinates for ic in initial_conditions)...)
-        velocity = hcat((ic.velocity for ic in initial_conditions)...)
-        mass = vcat((ic.mass for ic in initial_conditions)...)
-        density = vcat((ic.density for ic in initial_conditions)...)
-
-        # TODO: Throw warning when particles are overlapping
-        return new{eltype(coordinates)}(coordinates, velocity, mass, density)
+        return new{eltype(coordinates)}(particle_spacing, coordinates, velocities, masses,
+                                        densities)
     end
 end
 
@@ -42,3 +30,112 @@ end
 end
 
 @inline nparticles(initial_condition::InitialCondition) = length(initial_condition.mass)
+
+function Base.union(initial_condition::InitialCondition, initial_conditions...)
+    particle_spacing = initial_condition.particle_spacing
+    ic = first(initial_conditions)
+
+    if ndims(ic) != ndims(initial_condition)
+        throw(ArgumentError("all passed initial conditions must have the same dimensionality"))
+    end
+
+    if particle_spacing < eps()
+        throw(ArgumentError("all passed initial conditions must store a particle spacing"))
+    end
+
+    if !isapprox(ic.particle_spacing, particle_spacing)
+        throw(ArgumentError("all passed initial conditions must have the same particle spacing"))
+    end
+
+    too_close = find_too_close_particles(ic.coordinates, initial_condition.coordinates,
+                                         0.75particle_spacing)
+    valid_particles = setdiff(eachparticle(ic), too_close)
+
+    coordinates = hcat(initial_condition.coordinates, ic.coordinates[:, valid_particles])
+    velocity = hcat(initial_condition.velocity, ic.velocity[:, valid_particles])
+    mass = vcat(initial_condition.mass, ic.mass[valid_particles])
+    density = vcat(initial_condition.density, ic.density[valid_particles])
+
+    result = InitialCondition(coordinates, velocity, mass, density,
+                              particle_spacing=particle_spacing)
+
+    return union(result, Base.tail(initial_conditions)...)
+end
+
+Base.union(initial_condition::InitialCondition) = initial_condition
+
+function Base.setdiff(initial_condition::InitialCondition, initial_conditions...)
+    ic = first(initial_conditions)
+
+    if ndims(ic) != ndims(initial_condition)
+        throw(ArgumentError("all passed initial conditions must have the same dimensionality"))
+    end
+
+    particle_spacing = initial_condition.particle_spacing
+    if particle_spacing < eps()
+        throw(ArgumentError("the initial condition in the first argument must store a particle spacing"))
+    end
+
+    too_close = find_too_close_particles(initial_condition.coordinates, ic.coordinates,
+                                         0.75particle_spacing)
+    valid_particles = setdiff(eachparticle(initial_condition), too_close)
+
+    coordinates = initial_condition.coordinates[:, valid_particles]
+    velocity = initial_condition.velocity[:, valid_particles]
+    mass = initial_condition.mass[valid_particles]
+    density = initial_condition.density[valid_particles]
+
+    result = InitialCondition(coordinates, velocity, mass, density,
+                              particle_spacing=particle_spacing)
+
+    return setdiff(result, Base.tail(initial_conditions)...)
+end
+
+Base.setdiff(initial_condition::InitialCondition) = initial_condition
+
+function Base.intersect(initial_condition::InitialCondition, initial_conditions...)
+    ic = first(initial_conditions)
+
+    if ndims(ic) != ndims(initial_condition)
+        throw(ArgumentError("all passed initial conditions must have the same dimensionality"))
+    end
+
+    particle_spacing = initial_condition.particle_spacing
+    if particle_spacing < eps()
+        throw(ArgumentError("the initial condition in the first argument must store a particle spacing"))
+    end
+
+    too_close = find_too_close_particles(initial_condition.coordinates, ic.coordinates,
+                                         0.75particle_spacing)
+
+    coordinates = initial_condition.coordinates[:, too_close]
+    velocity = initial_condition.velocity[:, too_close]
+    mass = initial_condition.mass[too_close]
+    density = initial_condition.density[too_close]
+
+    result = InitialCondition(coordinates, velocity, mass, density,
+                              particle_spacing=particle_spacing)
+
+    return intersect(result, Base.tail(initial_conditions)...)
+end
+
+Base.intersect(initial_condition::InitialCondition) = initial_condition
+
+# Find particles in `coords1` that are closer than `max_distance` to any particle in `coords2`
+function find_too_close_particles(coords1, coords2, max_distance)
+    NDIMS = size(coords1, 1)
+    result = Int[]
+
+    nhs = SpatialHashingSearch{NDIMS}(max_distance, size(coords2, 2))
+    TrixiParticles.initialize!(nhs, coords2)
+
+    # We are modifying the vector `result`, so this cannot be parallel
+    TrixiParticles.for_particle_neighbor(coords1, coords2, nhs,
+                                         parallel=false) do particle, _, _, _
+        if !(particle in result)
+            append!(result, particle)
+        end
+    end
+
+    return result
+end
diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl
@@ -342,7 +342,21 @@ end
 # By default, loop over `each_moving_particle(system)`.
 @inline function for_particle_neighbor(f, system, neighbor_system,
                                        system_coords, neighbor_coords, neighborhood_search;
-                                       particles=each_moving_particle(system))
+                                       particles=each_moving_particle(system),
+                                       parallel=true)
+    for_particle_neighbor(f, system_coords, neighbor_coords, neighborhood_search,
+                          particles=particles, parallel=parallel)
+end
+
+@inline function for_particle_neighbor(f, system_coords, neighbor_coords,
+                                       neighborhood_search;
+                                       particles=axes(system_coords, 2), parallel=true)
+    for_particle_neighbor(f, system_coords, neighbor_coords, neighborhood_search, particles,
+                          Val(parallel))
+end
+
+@inline function for_particle_neighbor(f, system_coords, neighbor_coords,
+                                       neighborhood_search, particles, parallel::Val{true})
     @threaded for particle in particles
         for_particle_neighbor_inner(f, system_coords, neighbor_coords, neighborhood_search,
                                     particle)
@@ -351,6 +365,16 @@ end
     return nothing
 end
 
+@inline function for_particle_neighbor(f, system_coords, neighbor_coords,
+                                       neighborhood_search, particles, parallel::Val{false})
+    for particle in particles
+        for_particle_neighbor_inner(f, system_coords, neighbor_coords, neighborhood_search,
+                                    particle)
+    end
+
+    return nothing
+end
+
 # Use this function barrier and unpack inside to avoid passing closures to Polyester.jl
 # with @batch (@threaded).
 # Otherwise, @threaded does not work here with Julia ARM on macOS.

diff --git a/src/setups/rectangular_shape.jl b/src/setups/rectangular_shape.jl
@@ -61,7 +61,8 @@ function RectangularShape(particle_spacing, n_particles_per_dimension,
     densities = density * ones(ELTYPE, n_particles)
     masses = density * particle_spacing^NDIMS * ones(ELTYPE, n_particles)
 
-    return InitialCondition(coordinates, velocities, masses, densities)
+    return InitialCondition(coordinates, velocities, masses, densities,
+                            particle_spacing=particle_spacing)
 end
 
 function rectangular_shape_coords(particle_spacing, n_particles_per_dimension,

diff --git a/src/setups/rectangular_tank.jl b/src/setups/rectangular_tank.jl
@@ -97,13 +97,13 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real}
                                                                             particle_spacing,
                                                                             spacing_ratio)
 
-        boundary_coordinates, face_indices = initialize_boundaries(particle_spacing /
-                                                                   spacing_ratio,
+        boundary_spacing = particle_spacing / spacing_ratio
+        boundary_coordinates, face_indices = initialize_boundaries(boundary_spacing,
                                                                    tank_size_,
                                                                    n_boundaries_per_dim,
                                                                    n_layers, faces)
 
-        boundary_masses = boundary_density * (particle_spacing / spacing_ratio)^NDIMS *
+        boundary_masses = boundary_density * boundary_spacing^NDIMS *
                           ones(ELTYPE, size(boundary_coordinates, 2))
         boundary_densities = boundary_density * ones(ELTYPE, size(boundary_coordinates, 2))
         boundary_velocities = zeros(ELTYPE, size(boundary_coordinates))
@@ -116,7 +116,8 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real}
                                  zeros(NDIMS), fluid_density, init_velocity=init_velocity)
 
         boundary = InitialCondition(boundary_coordinates, boundary_velocities,
-                                    boundary_masses, boundary_densities)
+                                    boundary_masses, boundary_densities,
+                                    particle_spacing=boundary_spacing)
 
         return new{NDIMS, 2 * NDIMS, ELTYPE}(fluid, boundary, fluid_size_, tank_size_,
                                              faces, face_indices,

diff --git a/src/setups/sphere_shape.jl b/src/setups/sphere_shape.jl
@@ -112,7 +112,8 @@ function SphereShape(particle_spacing, radius, center_position, density;
     masses = density * particle_spacing^NDIMS * ones(ELTYPE, n_particles)
     velocities = init_velocity .* ones(ELTYPE, size(coordinates))
 
-    return InitialCondition(coordinates, velocities, masses, densities)
+    return InitialCondition(coordinates, velocities, masses, densities,
+                            particle_spacing=particle_spacing)
 end
 
 """