Merge branch 'main' into dev

svchb · web-flow · commit fc62e47f6127 · 2025-08-06T15:00:48.000+02:00
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -2,6 +2,7 @@
 Asciicast = "2600d445-abca-43b9-92aa-ce144ac0b05b"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
diff --git a/docs/literate/src/tut_packing.jl b/docs/literate/src/tut_packing.jl
diff --git a/docs/make.jl b/docs/make.jl
@@ -3,6 +3,7 @@ using TrixiParticles
 using TrixiParticles.TrixiBase
 using TrixiParticles.PointNeighbors
 using Asciicast: Asciicast
+using Literate: Literate
 
 # Get TrixiParticles.jl root directory
 trixiparticles_root_dir = dirname(@__DIR__)
@@ -83,6 +84,9 @@ replace_with_code(joinpath("docs", "src", "tutorials_template", "tut_dam_break.m
 replace_with_code(joinpath("docs", "src", "tutorials_template", "tut_beam.md"))
 replace_with_code(joinpath("docs", "src", "tutorials_template", "tut_falling.md"))
 
+Literate.markdown(joinpath("docs", "literate", "src", "tut_packing.jl"),
+                  joinpath("docs", "src", "tutorials"))
+
 copy_file("AUTHORS.md",
           "in the [LICENSE.md](LICENSE.md) file" => "under [License](@ref)")
 copy_file("CONTRIBUTING.md",
diff --git a/docs/src/preprocessing/preprocessing.md b/docs/src/preprocessing/preprocessing.md
@@ -302,7 +302,8 @@ Pages = [joinpath("preprocessing", "geometries", "io.jl")]
 ```
 
 # [Particle Packing](@id particle_packing)
-To obtain a body-fitted and isotropic particle distribution, an initial configuration (see [Sampling of Geometries](@ref sampling_of_geometries)) is first generated. This configuration is then packed using a [`ParticlePackingSystem`](@ref).
+To obtain a body-fitted and isotropic particle distribution, an initial configuration (see [Sampling of Geometries](@ref sampling_of_geometries)) is first generated. This configuration is then packed using a [`ParticlePackingSystem`](@ref) following the steps introduced in [neher2025robustefficientpreprocessingtechniques](@cite).
+For a hands-on tutorial with complete examples, see the particle packing tutorial.
 The preprocessing pipeline consists of the following steps:
 
 - Load geometry: Fig. 1, [`load_geometry`](@ref).
diff --git a/docs/src/refs.bib b/docs/src/refs.bib
@@ -587,6 +587,16 @@ @Article{Negi2022
   publisher = {AIP Publishing},
 }
 
+@misc{neher2025robustefficientpreprocessingtechniques,
+      title={Robust and efficient pre-processing techniques for particle-based methods including dynamic boundary generation},
+      author={Niklas S. Neher and Erik Faulhaber and Sven Berger and Christian Weißenfels and Gregor J. Gassner and Michael Schlottke-Lakemper},
+      year={2025},
+      eprint={2506.21206},
+      archivePrefix={arXiv},
+      primaryClass={math.NA},
+      url={https://arxiv.org/abs/2506.21206},
+}
+
 @Article{O’Connor2021,
   author    = {O’Connor, Joseph and Rogers, Benedict D.},
   journal   = {Journal of Fluids and Structures},
diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md
@@ -1,4 +1,4 @@
-# Tutorial
+# Tutorials
 
 
 ## Fluid
@@ -11,3 +11,6 @@
 
 ## Fluid-Structure Interaction
 - [Setting up a falling structure](tutorials/tut_falling.md)
+
+## Preprocessing
+- [Particle packing](tutorials/tut_packing.md)
diff --git a/examples/preprocessing/data/potato.asc b/examples/preprocessing/data/potato.asc
@@ -0,0 +1,15 @@
+# ASCII
+1.14991 -0.00561128 0
+0.512738 -0.0262819 0
+0.18511 -0.134802 0
+0.00682607 -0.357011 0
+0.126199 -0.549764 0
+0.404735 -0.687223 0
+0.805744 -0.699626 0
+1.26308 -0.934754 0
+1.60518 -0.940955 0
+1.94779 -0.809697 0
+1.9943 -0.532194 0
+1.8031 -0.271745 0
+1.56642 -0.131702 0
+1.14991 -0.00561128 0
diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
@@ -56,8 +56,8 @@ include("callbacks/callbacks.jl")
 include("general/semidiscretization.jl")
 include("general/gpu.jl")
 include("io/io.jl")
-include("visualization/recipes_plots.jl")
 include("preprocessing/preprocessing.jl")
+include("visualization/recipes_plots.jl")
 
 export Semidiscretization, semidiscretize, restart_with!
 export InitialCondition
diff --git a/src/callbacks/steady_state_reached.jl b/src/callbacks/steady_state_reached.jl
@@ -26,17 +26,19 @@ end
 
 function SteadyStateReachedCallback(; interval::Integer=0, dt=0.0,
                                     interval_size::Integer=10, abstol=1.0e-8, reltol=1.0e-6)
+    ELTYPE = eltype(abstol)
     abstol, reltol = promote(abstol, reltol)
 
     if dt > 0 && interval > 0
         throw(ArgumentError("setting both `interval` and `dt` is not supported"))
     end
 
     if dt > 0
-        interval = Float64(dt)
+        interval = convert(ELTYPE, dt)
     end
 
-    steady_state_callback = SteadyStateReachedCallback(interval, abstol, reltol, [Inf64],
+    steady_state_callback = SteadyStateReachedCallback(interval, abstol, reltol,
+                                                       [convert(ELTYPE, Inf)],
                                                        interval_size)
 
     if dt > 0
diff --git a/src/general/custom_quantities.jl b/src/general/custom_quantities.jl
@@ -6,22 +6,26 @@ Returns the total kinetic energy of all particles in a system.
 function kinetic_energy(system, v_ode, u_ode, semi, t)
     v = wrap_v(v_ode, system, semi)
 
-    # If `each_moving_particle` is empty (no moving particles), return zero
-    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)
+    # TODO: `current_velocity` should only contain active particles
+    # (see https://github.com/trixi-framework/TrixiParticles.jl/issues/850)
+    velocity = reinterpret(reshape, SVector{ndims(system), eltype(v)},
+                           view(current_velocity(v, system), :, active_particles(system)))
+    mass = view(system.mass, active_particles(system))
+
+    return mapreduce(+, velocity, mass) do v_i, m_i
+        return m_i * dot(v_i, v_i) / 2
     end
 end
 
+kinetic_energy(system::BoundarySystem, v_ode, u_ode, semi, t) = zero(eltype(system))
+
 """
     total_mass
 
 Returns the total mass of all particles in a system.
 """
 function total_mass(system, v_ode, u_ode, semi, t)
-    return sum(eachparticle(system)) do particle
-        return system.mass[particle]
-    end
+    return sum(system.mass)
 end
 
 function total_mass(system::BoundarySystem, v_ode, u_ode, semi, t)
diff --git a/src/setups/complex_shape.jl b/src/setups/complex_shape.jl
@@ -78,8 +78,8 @@ function ComplexShape(geometry; particle_spacing, density,
 end
 
 """
-    sample_boundary(signed_distance_field::SignedDistanceField;
-                    boundary_thickness::Real, place_on_shell=true)
+    sample_boundary(signed_distance_field;
+                    boundary_density, boundary_thickness, tlsph=true)
 
 Sample boundary particles of a complex geometry by using the [`SignedDistanceField`](@ref)
 of the geometry.
@@ -90,6 +90,7 @@ of the geometry.
 # Keywords
 - `boundary_thickness`: Thickness of the boundary
 - `place_on_shell`:     When `place_on_shell=true`, boundary particles will be placed
+- `boundary_density`: Density of each boundary particle.
                         one particle spacing from the surface of the geometry.
                         Otherwise when `place_on_shell=true` (simulating fluid particles),
                         boundary particles will be placed half particle spacing away from the surface.
diff --git a/src/visualization/recipes_plots.jl b/src/visualization/recipes_plots.jl
@@ -147,3 +147,26 @@ RecipesBase.@recipe function f(::Union{InitialCondition, Semidiscretization},
         end
     end
 end
+
+RecipesBase.@recipe function f(geometry::Polygon)
+    x = stack(geometry.vertices)[1, :]
+    y = stack(geometry.vertices)[2, :]
+
+    aspect_ratio --> :equal
+    grid --> false
+
+    # First plot the vertices as a scatter plot
+    @series begin
+        seriestype --> :scatter
+        return (x, y)
+    end
+
+    # Now plot the edges as a line plot
+    @series begin
+        # Ignore series in legend and color cycling. Note that `:=` forces the attribute,
+        # whereas `-->` would only set it if it is not already set.
+        primary := false
+
+        return (x, y)
+    end
+end
diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl
@@ -354,6 +354,44 @@ end
             backend = TrixiParticles.KernelAbstractions.get_backend(sol.u[end].x[1])
             @test backend == Main.parallelization_backend
         end
+
+        @trixi_testset "fluid/pipe_flow_2d.jl - steady state reached (`dt`)" begin
+            steady_state_reached = SteadyStateReachedCallback(; dt=0.002f0,
+                                                              interval_size=10,
+                                                              reltol=1.0f-3)
+
+            @trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__,
+                                                             joinpath(examples_dir(),
+                                                                      "fluid",
+                                                                      "pipe_flow_2d.jl"),
+                                                             extra_callback=steady_state_reached,
+                                                             tspan=(0.0f0, 1.5f0),
+                                                             parallelization_backend=Main.parallelization_backend,
+                                                             viscosity_boundary=nothing)
+
+            # Make sure that the simulation is terminated after a reasonable amount of time
+            @test 0.1 < sol.t[end] < 1.0
+            @test sol.retcode == ReturnCode.Terminated
+        end
+
+        @trixi_testset "fluid/pipe_flow_2d.jl - steady state reached (`interval`)" begin
+            steady_state_reached = SteadyStateReachedCallback(; interval=1,
+                                                              interval_size=10,
+                                                              reltol=1.0f-3)
+            @trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__,
+                                                             joinpath(examples_dir(),
+                                                                      "fluid",
+                                                                      "pipe_flow_2d.jl"),
+                                                             extra_callback=steady_state_reached,
+                                                             dtmax=2.0f-3,
+                                                             tspan=(0.0f0, 1.5f0),
+                                                             parallelization_backend=Main.parallelization_backend,
+                                                             viscosity_boundary=nothing)
+
+            # Make sure that the simulation is terminated after a reasonable amount of time
+            @test 0.1 < sol.t[end] < 1.0
+            @test sol.retcode == ReturnCode.Terminated
+        end
     end
 
     @testset verbose=true "Solid" begin
diff --git a/test/general/custom_quantities.jl b/test/general/custom_quantities.jl
@@ -0,0 +1,115 @@
+@testset verbose=true "Custom Quantities" begin
+    particle_spacing = 0.1
+    coordinates = [0.0 0.1 0.2; 0.0 0.0 0.1]
+    velocities = [1.0 2.0 0.5; 0.5 1.0 1.5]
+    densities = [1000.0, 1000.0, 1000.0]
+    pressures = [101325.0, 101330.0, 101320.0]
+
+    initial_condition = InitialCondition(; coordinates, velocity=velocities,
+                                         density=densities,
+                                         pressure=pressures, particle_spacing)
+
+    smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+    smoothing_length = 1.2 * particle_spacing
+
+    fluid_system = EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel,
+                                               smoothing_length, 1.0)
+    fluid_system.cache.density .= initial_condition.density
+
+    boundary_model = BoundaryModelDummyParticles(initial_condition.density,
+                                                 initial_condition.mass,
+                                                 AdamiPressureExtrapolation(),
+                                                 smoothing_kernel, smoothing_length)
+
+    boundary_system = BoundarySPHSystem(initial_condition, boundary_model)
+
+    semi = Semidiscretization(fluid_system, boundary_system)
+
+    v_ode, u_ode = semidiscretize(semi, (0, 1)).u0.x
+    t = 0.0
+
+    @testset "Kinetic Energy" begin
+        @testset "Fluid System" begin
+            ekin = kinetic_energy(fluid_system, v_ode, u_ode, semi, t)
+
+            expected_ekin = sum(velocities) do velocity
+                return 0.5 * first(fluid_system.mass) * dot(velocity, velocity)
+            end
+
+            @test isapprox(expected_ekin, ekin)
+        end
+
+        @testset "Boundary System" begin
+            ekin = kinetic_energy(boundary_system, v_ode, u_ode, semi, t)
+            @test ekin == 0
+        end
+    end
+
+    @testset "Total Mass" begin
+        @testset "Fluid System" begin
+            mass = total_mass(fluid_system, v_ode, u_ode, semi, t)
+            expected_mass = first(fluid_system.mass) * nparticles(fluid_system)
+
+            @test isapprox(mass, expected_mass)
+        end
+
+        @testset "Boundary System" begin
+            mass = total_mass(boundary_system, v_ode, u_ode, semi, t)
+
+            @test isnan(mass)
+        end
+    end
+
+    @testset "Pressure Quantities" begin
+        @testset "Max Pressure" begin
+            max_p = max_pressure(fluid_system, v_ode, u_ode, semi, t)
+            @test isapprox(max_p, 101330.0)
+
+            # Boundary system should return NaN
+            @test isnan(max_pressure(boundary_system, v_ode, u_ode, semi, t))
+        end
+
+        @testset "Min Pressure" begin
+            min_p = min_pressure(fluid_system, v_ode, u_ode, semi, t)
+            @test min_p ≈ 101320.0
+
+            # Boundary system should return NaN
+            @test isnan(min_pressure(boundary_system, v_ode, u_ode, semi, t))
+        end
+
+        @testset "Average Pressure" begin
+            avg_p = avg_pressure(fluid_system, v_ode, u_ode, semi, t)
+            expected_avg = (101325.0 + 101330.0 + 101320.0) / 3
+            @test isapprox(avg_p, expected_avg)
+
+            # Boundary system should return NaN
+            @test isnan(avg_pressure(boundary_system, v_ode, u_ode, semi, t))
+        end
+    end
+
+    @testset "Density Quantities" begin
+        @testset "max_density" begin
+            max_d = max_density(fluid_system, v_ode, u_ode, semi, t)
+            @test isapprox(max_d, 1000.0)  # All particles have same density
+
+            # Boundary system should return NaN
+            @test isnan(max_density(boundary_system, v_ode, u_ode, semi, t))
+        end
+
+        @testset "min_density" begin
+            min_d = min_density(fluid_system, v_ode, u_ode, semi, t)
+            @test isapprox(min_d, 1000.0)
+
+            # Boundary system should return NaN
+            @test isnan(min_density(boundary_system, v_ode, u_ode, semi, t))
+        end
+
+        @testset "avg_density" begin
+            avg_d = avg_density(fluid_system, v_ode, u_ode, semi, t)
+            @test isapprox(avg_d, 1000.0)
+
+            # Boundary system should return NaN
+            @test isnan(avg_density(boundary_system, v_ode, u_ode, semi, t))
+        end
+    end
+end
diff --git a/test/general/general.jl b/test/general/general.jl
@@ -5,3 +5,4 @@ include("semidiscretization.jl")
 include("interpolation.jl")
 include("buffer.jl")
 include("util.jl")
+include("custom_quantities.jl")
