Particle packing (#529)

* visualize signed distances * fix infinite normals * clean up TriangleMesh * clean up TriangleMesh * modify shapes * change ELTYPE in IO * add files and remove `scale_factor` * add shape and `complexShape` tests * remove complex tank example * add docs for `complexShapes` * add docs for the algorithms * add test data * add tests for sampled shapes * data path relative to module * use `Threads.@threads` instead of `threaded` * fix @threaded * fix tests * fix wrong directory * use `@threaded` again * change winding number factor for julia 1.9 * fix tests * modify examples * revise * add docs for `ParticlePackingSystem` * add signed distance docs * add closure for sdf * modify example 3d * add setter for tvf * add trixi2vtk(sdf) * sample boundary in `ComplexShape` * add stl files * add `ComplexShape` structure * fix callback used check * implement suggestions * rename kwargs * add comments * apply formatter * implement suggestion * fix tests * remove abstract type `Shapes` * fix tests * make type stable * implement suggestions * add `unique_soted` * add comment * implement suggestions * RADME and NEWS entry * adapt docs * `for` instead of `while` * rename `shape` to `geometry` * rename again * implement suggestions * fix typo * implement suggestions * adapt example files * sum-init with `zero(eltype))` * implement suggestions * construct hierarchy in constructor * use ASCII characters * minor changes * add 2d hierarchical winding * minor changes * change function name * fix * implement suggestions * fix typo * apply formatter * fix typo * add `deleteat!` * circle example eps() -> 0 * start docs * simplify * implement suggestions * fix docs * fix tests * fix tests * minor changes * typo * fix * fix again * fx * fix tests * fix again * implement suggestions * modify NEWS.md * modify again * remove face nhs * adapt examples * minor changes * clean up * add tgv validation example * add ldc validation * modify validation * increase `maxiters` * remove double velocity initialization * add random displacement * implement suggestions * update `NEWS.md` * remove check for viscosity * fix tests * add short tilde description in docs * remove discarded example from tests * implement suggestions * apply formatter * add boundary compress factor * fix signed distance field * fix io * implement suggestions * add comment * fix test * fix typo * apply formatter * fix unit tests * apply formatter * revise docs * add autodocs * fix tests * implement suggestions * remove dependency * implement suggestions nhs * add `*` * formatting docs * improve `triangle_plane_intersection` * minor changes * implement sdf suggestions * fix kwarg bug * fix sign bit bug * add comment * revise docs * improve API * adapt example files * terminate with steady state callback * fix tests * implement suggestions * add unit test for `SignedDistanceField` * add system tests * apply formatter * remove geometry from system * add `delete_positions_in_empty_cells` * revise face nhs * fix face nhs * add cell intersection tests * fix typos * fix tests * add kwargs for trixi2vtk * resolve dependency error * implement suggestions * implement suggestions * add `neighbors` field in nhs * fix * add example * remove cell intersection check * adapt nhs tests * fix tests * fix doc tests * add comment for update * add whole packing dir to docs * implement suggestions * add comments * fix tests * add comments * fix rounding error for v1.9 * Revert "fix rounding error for v1.9" This reverts commit f026ddc. * implement suggestions * Increase Version * Update NEWS.md * Typo * Typo * Update NEWS.md Co-authored-by: Erik Faulhaber <[email protected]> * modify news --------- Co-authored-by: LasNikas <[email protected]> Co-authored-by: Sven Berger <[email protected]> Co-authored-by: Erik Faulhaber <[email protected]>
trixi-framework · Jan 17, 2025 · c36ac19 · svchb · Jan 17, 2025 · c36ac19
1 parent 7c537a3
commit c36ac19
Show file tree

Hide file tree

Showing 33 changed files with 7,576 additions and 85 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -3,51 +3,86 @@
 TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 
+## Version 0.2.5
+
+### Features
+- Add particle packing for 2D (.asc) and 3D (.stl) geometries (#529)
+
+### Compatibility Changes
+- Dropped support for Julia 1.9
+
+## Version 0.2.4
+
+### Features
+
+- A method to prevent penetration of fast moving particles with solids was added (#498)
+- Added the callback `SteadyStateReachedCallback` to detect convergence of static simulations (#601)
+- Added ideal gas state equation (#607)
+- Simulations can be run with `Float32` (#662)
+
+### Documentation
+
+- Documentation for GPU support was added (#660)
+- A new user tutorial was added (#514)
+- Several docs fixes (#663, #659, #637, #658, #664)
+
+
 ## Version 0.2.3
 
 ### Highlights
+
 Transport Velocity Formulation (TVF) based on the work of Ramachandran et al. "Entropically damped artiﬁcial compressibility for SPH" (2019) was added.
 
 ## Version 0.2.2
 
 ### Highlights
+
 Hotfix for threaded sampling of complex geometries.
 
 ## Version 0.2.1
 
 ### Highlights
+
 Particle sampling of complex geometries from `.stl` and `.asc` files.
 
 ## Version 0.2.0
 
 ### Removed
+
 Use of the internal neighborhood search has been removed and replaced with PointNeighbors.jl.
 
 ## Development Cycle 0.1
 
 ### Highlights
 
 #### Discrete Element Method
+
 A basic implementation of the discrete element method was added.
 
 #### Surface Tension and Adhesion Model
+
 A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids" (2013) was added to WCSPH.
 
 #### Support for Open Boundaries
+
 Open boundaries using the method of characteristics based on the work of Lastiwka et al., "Permeable and non-reflecting boundary conditions in SPH" (2009) were added for WCSPH and EDAC.
 
 ## Pre Initial Release (v0.1.0)
+
 This section summarizes the initial features that TrixiParticles.jl was released with.
 
 ### Highlights
 #### EDAC
+
 An implementation of EDAC (Entropically Damped Artificial Compressibility) was added,
 which allows for more stable simulations compared to basic WCSPH and reduces spurious pressure oscillations.
 
 #### WCSPH
+
 An implementation of WCSPH (Weakly Compressible Smoothed Particle Hydrodynamics), which is the classical SPH approach.
 
 Features:
+
 - Correction schemes (Shepard (0. Order) ... MixedKernelGradient (1. Order))
 - Density reinitialization
 - Kernel summation and Continuity equation density formulations
@@ -57,4 +92,5 @@ Features:
 
 
 #### TLSPH
+
 An implementation of TLSPH (Total Lagrangian Smoothed Particle Hydrodynamics) for solid bodies enabling FSI (Fluid Structure Interactions).
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 authors = ["erik.faulhaber <[email protected]>"]
-version = "0.2.5-dev"
+version = "0.2.5"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

diff --git a/docs/src/preprocessing/preprocessing.md b/docs/src/preprocessing/preprocessing.md
@@ -255,3 +255,12 @@ Pages = [joinpath("preprocessing", "point_in_poly", "winding_number_jacobson.jl"
 Modules = [TrixiParticles]
 Pages = [joinpath("preprocessing", "geometries", "io.jl")]
 ```
+
+# Particle Packing
+
+TODO
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = map(file -> joinpath("preprocessing", "particle_packing", file), readdir(joinpath("..", "src", "preprocessing", "particle_packing")))
+```
diff --git a/examples/preprocessing/complex_shape_2d.jl b/examples/preprocessing/complex_shape_2d.jl
@@ -21,6 +21,7 @@ shape_sampled = ComplexShape(geometry; particle_spacing, density=1.0,
 
 trixi2vtk(shape_sampled.initial_condition)
 
+# trixi2vtk(shape_sampled.signed_distance_field)
 # trixi2vtk(shape_sampled.grid, w=shape_sampled.winding_numbers)
 
 # Plot the winding number field

diff --git a/examples/preprocessing/complex_shape_3d.jl b/examples/preprocessing/complex_shape_3d.jl
@@ -13,6 +13,9 @@ point_in_geometry_algorithm = WindingNumberJacobson(; geometry,
 
 # Returns `InitialCondition`
 shape_sampled = ComplexShape(geometry; particle_spacing, density=1.0,
-                             point_in_geometry_algorithm)
+                             boundary_thickness=5 * particle_spacing,
+                             create_signed_distance_field=true,
+                             sample_boundary=false, point_in_geometry_algorithm)
 
-trixi2vtk(shape_sampled)
+trixi2vtk(shape_sampled.initial_condition, filename="initial_condition_fluid")
+trixi2vtk(shape_sampled.signed_distance_field)
diff --git a/examples/preprocessing/packing_2d.jl b/examples/preprocessing/packing_2d.jl
@@ -0,0 +1,89 @@
+using TrixiParticles
+using OrdinaryDiffEq
+
+filename = "circle"
+file = pkgdir(TrixiParticles, "examples", "preprocessing", "data", filename * ".asc")
+
+# ==========================================================================================
+# ==== Packing parameters
+save_intervals = false
+tlsph = true
+
+# ==========================================================================================
+# ==== Resolution
+particle_spacing = 0.03
+
+# The following depends on the sampling of the particles. In this case `boundary_thickness`
+# means literally the thickness of the boundary packed with boundary particles and *not*
+# how many rows of boundary particles will be sampled.
+boundary_thickness = 8particle_spacing
+
+# ==========================================================================================
+# ==== Load complex geometry
+density = 1000.0
+
+geometry = load_geometry(file)
+
+signed_distance_field = SignedDistanceField(geometry, particle_spacing;
+                                            use_for_boundary_packing=true,
+                                            max_signed_distance=boundary_thickness)
+
+point_in_geometry_algorithm = WindingNumberJacobson(; geometry,
+                                                    winding_number_factor=0.4,
+                                                    hierarchical_winding=true)
+# Returns `InitialCondition`
+shape_sampled = ComplexShape(geometry; particle_spacing, density,
+                             point_in_geometry_algorithm)
+
+# Returns `InitialCondition`
+boundary_sampled = sample_boundary(signed_distance_field; boundary_density=density,
+                                   boundary_thickness, tlsph)
+
+trixi2vtk(shape_sampled)
+trixi2vtk(boundary_sampled, filename="boundary")
+
+# ==========================================================================================
+# ==== Packing
+
+# Large `background_pressure` can cause high accelerations. That is, the adaptive
+# time-stepsize will be adjusted properly. We found that the following order of
+# `background_pressure` result in appropriate stepsizes.
+background_pressure = 1e6 * particle_spacing^ndims(geometry)
+
+packing_system = ParticlePackingSystem(shape_sampled;
+                                       signed_distance_field, tlsph=tlsph,
+                                       background_pressure)
+
+boundary_system = ParticlePackingSystem(boundary_sampled;
+                                        is_boundary=true, signed_distance_field,
+                                        tlsph=tlsph, boundary_compress_factor=0.8,
+                                        background_pressure)
+
+# ==========================================================================================
+# ==== Simulation
+semi = Semidiscretization(packing_system, boundary_system)
+
+# Use a high `tspan` to guarantee that the simulation runs at least for `maxiters`
+tspan = (0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+# Use this callback to stop the simulation when it is sufficiently close to a steady state
+steady_state = SteadyStateReachedCallback(; interval=1, interval_size=10,
+                                          abstol=1.0e-5, reltol=1.0e-3)
+
+info_callback = InfoCallback(interval=50)
+
+saving_callback = save_intervals ?
+                  SolutionSavingCallback(interval=10, prefix="", ekin=kinetic_energy) :
+                  nothing
+
+callbacks = CallbackSet(UpdateCallback(), saving_callback, info_callback, steady_state)
+
+sol = solve(ode, RDPK3SpFSAL35();
+            save_everystep=false, maxiters=1000, callback=callbacks, dtmax=1e-2)
+
+packed_ic = InitialCondition(sol, packing_system, semi)
+packed_boundary_ic = InitialCondition(sol, boundary_system, semi)
+
+trixi2vtk(packed_ic, filename="initial_condition_packed")
+trixi2vtk(packed_boundary_ic, filename="initial_condition_boundary_packed")
diff --git a/examples/preprocessing/packing_3d.jl b/examples/preprocessing/packing_3d.jl
@@ -0,0 +1,24 @@
+using TrixiParticles
+using OrdinaryDiffEq
+
+filename = "sphere"
+file = pkgdir(TrixiParticles, "examples", "preprocessing", "data", filename * ".stl")
+
+# ==========================================================================================
+# ==== Packing parameters
+tlsph = true
+save_intervals = false
+
+# ==========================================================================================
+# ==== Resolution
+particle_spacing = 0.1
+
+# The following depends on the sampling of the particles. In this case `boundary_thickness`
+# means literally the thickness of the boundary packed with boundary particles and *not*
+# how many rows of boundary particles will be sampled.
+boundary_thickness = 8particle_spacing
+
+trixi_include(joinpath(examples_dir(), "preprocessing", "packing_2d.jl"),
+              density=1000.0, particle_spacing=particle_spacing, file=file,
+              boundary_thickness=boundary_thickness, tlsph=tlsph,
+              save_intervals=save_intervals)
diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
@@ -40,7 +40,6 @@ using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, paraview_collection, vtk_save
 include("general/system.jl")
 # `util.jl` needs to be next because of the macros `@trixi_timeit` and `@threaded`
 include("util.jl")
-include("preprocessing/preprocessing.jl")
 include("callbacks/callbacks.jl")
 include("general/general.jl")
 include("setups/setups.jl")
@@ -52,6 +51,7 @@ include("general/semidiscretization.jl")
 include("general/gpu.jl")
 include("visualization/write2vtk.jl")
 include("visualization/recipes_plots.jl")
+include("preprocessing/preprocessing.jl")
 
 export Semidiscretization, semidiscretize, restart_with!
 export InitialCondition
@@ -77,8 +77,10 @@ export BoundaryMovement
 export examples_dir, validation_dir, trixi_include
 export trixi2vtk
 export RectangularTank, RectangularShape, SphereShape, ComplexShape
+export ParticlePackingSystem, SignedDistanceField
 export WindingNumberHormann, WindingNumberJacobson
-export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry
+export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry,
+       sample_boundary
 export SourceTermDamping
 export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection,
        GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection

diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl
@@ -85,6 +85,10 @@ function (update_callback!::UpdateCallback)(integrator)
         update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t)
     end
 
+    @trixi_timeit timer() "update particle packing" foreach_system(semi) do system
+        update_particle_packing(system, v_ode, u_ode, semi, integrator)
+    end
+
     @trixi_timeit timer() "update TVF" foreach_system(semi) do system
         update_transport_velocity!(system, v_ode, semi)
     end

diff --git a/src/general/custom_quantities.jl b/src/general/custom_quantities.jl
@@ -5,7 +5,7 @@ Returns the total kinetic energy of all particles in a system.
 """
 function kinetic_energy(v, u, t, system)
     # If `each_moving_particle` is empty (no moving particles), return zero
-    return sum(each_moving_particle(system), init=0.0) do particle
+    return sum(each_moving_particle(system), init=zero(eltype(system))) do particle
         velocity = current_velocity(v, system, particle)
         return 0.5 * system.mass[particle] * dot(velocity, velocity)
     end

diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl
@@ -326,6 +326,37 @@ end
 
 Base.intersect(initial_condition::InitialCondition) = initial_condition
 
+function InitialCondition(sol::ODESolution, system, semi; use_final_velocity=false,
+                          min_particle_distance=system.initial_condition.particle_spacing /
+                                                4)
+    ic = system.initial_condition
+
+    v_ode, u_ode = sol.u[end].x
+
+    u = wrap_u(u_ode, system, semi)
+    v = wrap_u(v_ode, system, semi)
+
+    # Check if particles come too close especially when the surface exhibits large curvature
+    too_close = find_too_close_particles(u, min_particle_distance)
+
+    velocity_ = use_final_velocity ? view(v, 1:ndims(system), :) : ic.velocity
+
+    not_too_close = setdiff(eachparticle(system), too_close)
+
+    coordinates = u[:, not_too_close]
+    velocity = velocity_[:, not_too_close]
+    mass = ic.mass[not_too_close]
+    density = ic.density[not_too_close]
+    pressure = ic.pressure[not_too_close]
+
+    if length(too_close) > 0
+        @info "Removed $(length(too_close)) particles that are too close together"
+    end
+
+    return InitialCondition{ndims(ic)}(coordinates, velocity, mass, density, pressure,
+                                       ic.particle_spacing)
+end
+
 # 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)
@@ -338,7 +369,27 @@ function find_too_close_particles(coords1, coords2, max_distance)
     # We are modifying the vector `result`, so this cannot be parallel
     foreach_point_neighbor(coords1, coords2, nhs, parallel=false) do particle, _, _, _
         if !(particle in result)
-            append!(result, particle)
+            push!(result, particle)
+        end
+    end
+
+    return result
+end
+
+# Find particles in `coords` that are closer than `min_distance` to any other particle in `coords`
+function find_too_close_particles(coords, min_distance)
+    NDIMS = size(coords, 1)
+    result = Int[]
+
+    nhs = GridNeighborhoodSearch{NDIMS}(search_radius=min_distance,
+                                        n_points=size(coords, 2))
+    TrixiParticles.initialize!(nhs, coords)
+
+    # We are modifying the vector `result`, so this cannot be parallel
+    foreach_point_neighbor(coords, coords, nhs, parallel=false) do particle, neighbor, _, _
+        # Only consider particles with neighbors that are not to be removed
+        if particle != neighbor && !(particle in result) && !(neighbor in result)
+            push!(result, particle)
         end
     end