Refactor SPH system initialization to use keyword arguments

README.md

@@ -133,7 +133,7 @@ You can find the documentation for the latest release
   [![arXiv:2506.21206](https://img.shields.io/badge/arXiv-2506.21206-yellow)](https://arxiv.org/abs/2506.21206)
   [![doi:10.1016/j.cpc.2025.109898](https://zenodo.org/badge/doi/10.1016/j.cpc.2025.109898.svg)](https://doi.org/10.1016/j.cpc.2025.109898)
   [![reproduce me!](https://img.shields.io/badge/reproduce-me!-brightgreen)](https://github.com/trixi-framework/paper-2025-particle-based_preprocessing)
-  
+
 ## Cite Us
 
 If you use TrixiParticles.jl in your own research or write a paper using results obtained

docs/literate/src/tut_custom_kernel.jl

@@ -117,13 +117,12 @@ nothing # hide
 # `examples/fluid/dam_break_2d.jl` with our custom kernel and the corresponding
 # smoothing length.
 # ```@cast @__NAME__; width=100, height=50, delay=0, loop=true, loop_delay=5
-# trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
+# trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl");
 #               smoothing_kernel=MyGaussianKernel(),
-#               smoothing_length=smoothing_length);
+#               smoothing_length);
 # ```
-trixi_include(joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), #!md
-              smoothing_kernel=MyGaussianKernel(), #!md
-              smoothing_length=smoothing_length) #!md
+trixi_include(joinpath(examples_dir(), "fluid", "dam_break_2d.jl"); #!md
+              smoothing_kernel=MyGaussianKernel(), smoothing_length) #!md
 
 # See [Visualization](@ref) for how to visualize the final solution.
 # For the simplest visualization, we can use [Plots.jl](https://docs.juliaplots.org/stable/):

docs/literate/src/tut_packing.jl

@@ -48,8 +48,8 @@ signed_distance_field = SignedDistanceField(geometry, particle_spacing;
 
 # We can also visualize the SDF by simply creating an `InitialCondition` from the sampled points and plot it.
 # The color coding represents the signed distance to the geometry surface.
-sdf_ic = InitialCondition(; coordinates=stack(signed_distance_field.positions),
-                          density=1.0, particle_spacing=particle_spacing)
+sdf_ic = InitialCondition(; coordinates=stack(signed_distance_field.positions), density=1.0,
+                          particle_spacing)
 
 plot(sdf_ic, zcolor=signed_distance_field.distances, label=nothing, color=:coolwarm)
 plot!(geometry, linestyle=:dash, label=nothing, showaxis=false, color=:black,
@@ -63,8 +63,8 @@ signed_distance_field = SignedDistanceField(geometry, particle_spacing;
                                             max_signed_distance=boundary_thickness)
 
 # We can see in the plot that the SDF has been extended outwards to twice `max_signed_distance`.
-sdf_ic = InitialCondition(; coordinates=stack(signed_distance_field.positions),
-                          density=1.0, particle_spacing=particle_spacing)
+sdf_ic = InitialCondition(; coordinates=stack(signed_distance_field.positions), density=1.0,
+                          particle_spacing)
 
 plot(sdf_ic, zcolor=signed_distance_field.distances, label=nothing, color=:coolwarm)
 
@@ -105,7 +105,7 @@ point_in_geometry_algorithm = WindingNumberJacobson(; geometry)
 # This function creates an [`InitialCondition`](@ref InitialCondition) for the interior particles.
 # Here, we need to specify `particle_spacing` and `density`. For further arguments, please refer to
 # the documentation of [ComplexShape](@ref ComplexShape) or [InitialCondition](@ref InitialCondition).
-shape_sampled = ComplexShape(geometry; particle_spacing, density=density,
+shape_sampled = ComplexShape(geometry; particle_spacing, density,
                              point_in_geometry_algorithm)
 
 # If we want to assign the mass of each sampled particle consistently with its density,
@@ -146,9 +146,7 @@ smoothing_length = 0.8 * particle_spacing
 
 # Now we can create the packing system. For learning purposes, let’s first try
 # passing no signed distance field (SDF) and see what happens.
-packing_system = ParticlePackingSystem(shape_sampled;
-                                       smoothing_kernel=smoothing_kernel,
-                                       smoothing_length=smoothing_length,
+packing_system = ParticlePackingSystem(shape_sampled; smoothing_kernel, smoothing_length,
                                        signed_distance_field=nothing, background_pressure)
 
 # We now proceed with the familiar steps
@@ -163,8 +161,7 @@ maxiters = 100
 callbacks = CallbackSet(UpdateCallback())
 time_integrator = RDPK3SpFSAL35()
 
-sol = solve(ode, time_integrator;
-            abstol=1e-7, reltol=1e-4, save_everystep=false, maxiters=maxiters,
+sol = solve(ode, time_integrator; abstol=1e-7, reltol=1e-4, save_everystep=false, maxiters,
             callback=callbacks)
 
 packed_ic = InitialCondition(sol, packing_system, semi)
@@ -176,11 +173,8 @@ plot!(geometry, seriestype=:path, linewidth=2, color=:black, label=nothing)
 # geometric surface.
 
 # We therefore add an SDF for the geometry and repeat the same procedure.
-packing_system = ParticlePackingSystem(shape_sampled;
-                                       smoothing_kernel=smoothing_kernel,
-                                       smoothing_length=smoothing_length,
-                                       signed_distance_field,
-                                       background_pressure)
+packing_system = ParticlePackingSystem(shape_sampled; smoothing_kernel, smoothing_length,
+                                       signed_distance_field, background_pressure)
 
 # Again, we follow the same steps for semidiscretization and time integration.
 semi = Semidiscretization(packing_system)
@@ -192,8 +186,7 @@ maxiters = 1000
 callbacks = CallbackSet(UpdateCallback())
 time_integrator = RDPK3SpFSAL35()
 
-sol = solve(ode, time_integrator;
-            abstol=1e-7, reltol=1e-4, save_everystep=false, maxiters=maxiters,
+sol = solve(ode, time_integrator; abstol=1e-7, reltol=1e-4, save_everystep=false, maxiters,
             callback=callbacks)
 
 packed_ic = InitialCondition(sol, packing_system, semi)
@@ -214,12 +207,10 @@ plot!(geometry, seriestype=:path, color=:black, label=nothing, linewidth=2)
 # A `boundary_compress_factor` of `0.8` or `0.9` works well for most shapes.
 # Since we have a relatively large particle spacing compared to the
 # geometry size in this example, we will choose `0.7`.
-boundary_system = ParticlePackingSystem(boundary_sampled;
-                                        is_boundary=true,
-                                        smoothing_kernel=smoothing_kernel,
-                                        smoothing_length=smoothing_length,
-                                        boundary_compress_factor=0.7,
-                                        signed_distance_field, background_pressure)
+boundary_system = ParticlePackingSystem(boundary_sampled; is_boundary=true,
+                                        smoothing_kernel, smoothing_length,
+                                        boundary_compress_factor=0.7, signed_distance_field,
+                                        background_pressure)
 
 # We can now couple the boundary system with the interior system:
 semi = Semidiscretization(packing_system, boundary_system)
@@ -231,8 +222,7 @@ maxiters = 1000
 callbacks = CallbackSet(UpdateCallback())
 time_integrator = RDPK3SpFSAL35()
 
-sol = solve(ode, time_integrator;
-            abstol=1e-7, reltol=1e-4, save_everystep=false, maxiters=maxiters,
+sol = solve(ode, time_integrator; abstol=1e-7, reltol=1e-4, save_everystep=false, maxiters,
             callback=callbacks)
 
 packed_ic = InitialCondition(sol, packing_system, semi)
@@ -254,11 +244,8 @@ plot!(geometry, seriestype=:path, color=:black, linestyle=:dash, linewidth=2, la
 # Because we are satisfied with this particle distribution, we want it to remain unchanged.
 # Therefore, we set `fixed_system=true` so that this system is not integrated further
 # but instead serves as a static boundary for the packing of other domains.
-fixed_system = ParticlePackingSystem(packed_ic;
-                                     smoothing_kernel=smoothing_kernel,
-                                     smoothing_length=smoothing_length,
-                                     signed_distance_field=nothing,
-                                     background_pressure,
+fixed_system = ParticlePackingSystem(packed_ic; smoothing_kernel, smoothing_length,
+                                     signed_distance_field=nothing, background_pressure,
                                      fixed_system=true)
 
 # Now we define a rectangular domain that we want to pack.
@@ -273,10 +260,8 @@ sampled_outer_domain = setdiff(tank_domain.fluid, packed_ic)
 plot(sampled_outer_domain, packed_ic)
 
 # Next, we create a packing system for the outer domain.
-packing_system = ParticlePackingSystem(sampled_outer_domain;
-                                       smoothing_kernel=smoothing_kernel,
-                                       smoothing_length=smoothing_length,
-                                       signed_distance_field=nothing,
+packing_system = ParticlePackingSystem(sampled_outer_domain; smoothing_kernel,
+                                       smoothing_length, signed_distance_field=nothing,
                                        background_pressure)
 
 # Since we do not want to sample a boundary for the outer domain,
@@ -293,8 +278,8 @@ maxiters = 1000
 callbacks = CallbackSet(UpdateCallback())
 time_integrator = RDPK3SpFSAL35()
 
-sol_1 = solve(ode, time_integrator; abstol=1e-7, reltol=1e-4,
-              save_everystep=false, maxiters=maxiters, callback=callbacks)
+sol_1 = solve(ode, time_integrator; abstol=1e-7, reltol=1e-4, save_everystep=false,
+              maxiters, callback=callbacks)
 
 packed_outer_domain = InitialCondition(sol_1, packing_system, semi)
 
@@ -323,24 +308,17 @@ plot(pack_domain, fixed_domain, packed_ic)
 
 # We can now treat the particles outside the window, along with the already
 # finalized configuration of the complex geometry, as fixed systems:
-fixed_system_1 = ParticlePackingSystem(fixed_domain;
-                                       smoothing_kernel=smoothing_kernel,
-                                       smoothing_length=smoothing_length,
-                                       signed_distance_field=nothing,
-                                       background_pressure, fixed_system=true)
-
-fixed_system_2 = ParticlePackingSystem(packed_ic;
-                                       smoothing_kernel=smoothing_kernel,
-                                       smoothing_length=smoothing_length,
-                                       signed_distance_field=nothing,
-                                       background_pressure, fixed_system=true)
+fixed_system_1 = ParticlePackingSystem(fixed_domain; smoothing_kernel, smoothing_length,
+                                       signed_distance_field=nothing, background_pressure,
+                                       fixed_system=true)
+
+fixed_system_2 = ParticlePackingSystem(packed_ic; smoothing_kernel, smoothing_length,
+                                       signed_distance_field=nothing, background_pressure,
+                                       fixed_system=true)
 
 # The window that we want to pack is passed to a moving packing system:
-packing_system = ParticlePackingSystem(pack_domain;
-                                       smoothing_kernel=smoothing_kernel,
-                                       smoothing_length=smoothing_length,
-                                       signed_distance_field=nothing,
-                                       background_pressure)
+packing_system = ParticlePackingSystem(pack_domain; smoothing_kernel, smoothing_length,
+                                       signed_distance_field=nothing, background_pressure)
 
 semi = Semidiscretization(packing_system, fixed_system_1, fixed_system_2)
 
@@ -351,8 +329,8 @@ maxiters = 1000
 callbacks = CallbackSet(UpdateCallback())
 time_integrator = RDPK3SpFSAL35()
 
-sol_2 = solve(ode, time_integrator; abstol=1e-7, reltol=1e-4,
-              save_everystep=false, maxiters=maxiters, callback=callbacks)
+sol_2 = solve(ode, time_integrator; abstol=1e-7, reltol=1e-4, save_everystep=false,
+              maxiters, callback=callbacks)
 
 packed_fluid_domain = InitialCondition(sol_2, packing_system, semi)
 

docs/literate/src/tut_rigid_body_fsi.jl

@@ -42,10 +42,10 @@ sound_speed = 100.0
 state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
                                    exponent=1)
 
-tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density,
-                       n_layers=boundary_layers, spacing_ratio=spacing_ratio,
-                       faces=(true, true, true, false),
-                       acceleration=(0.0, -gravity), state_equation=state_equation)
+tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density;
+                       n_layers=boundary_layers, spacing_ratio,
+                       faces=(true, true, true, false), acceleration=(0.0, -gravity),
+                       state_equation)
 nothing # hide
 
 # ## Rigid body geometry
@@ -97,10 +97,11 @@ fluid_density_calculator = ContinuityDensity()
 viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
 density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
 
-fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
-                                           state_equation, fluid_smoothing_kernel,
-                                           fluid_smoothing_length, viscosity=viscosity,
-                                           density_diffusion=density_diffusion,
+fluid_system = WeaklyCompressibleSPHSystem(tank.fluid;
+                                           smoothing_kernel=fluid_smoothing_kernel,
+                                           smoothing_length=fluid_smoothing_length,
+                                           density_calculator=fluid_density_calculator,
+                                           state_equation, viscosity, density_diffusion,
                                            acceleration=(0.0, -gravity))
 nothing # hide
 
@@ -166,10 +167,10 @@ nothing # hide
 boundary_density_calculator = AdamiPressureExtrapolation()
 tank_boundary_model = BoundaryModelDummyParticles(tank.boundary.density,
                                                   tank.boundary.mass,
-                                                  state_equation=state_equation,
                                                   boundary_density_calculator,
                                                   fluid_smoothing_kernel,
-                                                  fluid_smoothing_length)
+                                                  fluid_smoothing_length;
+                                                  state_equation)
 
 boundary_system = WallBoundarySystem(tank.boundary, tank_boundary_model)
 nothing # hide
@@ -180,10 +181,10 @@ function rigid_body_boundary_model(shape)
                           structure_particle_spacing^ndims(fluid_system)
 
     return BoundaryModelDummyParticles(hydrodynamic_densities, hydrodynamic_masses,
-                                       state_equation=state_equation,
                                        boundary_density_calculator,
                                        fluid_smoothing_kernel,
-                                       fluid_smoothing_length)
+                                       fluid_smoothing_length;
+                                       state_equation)
 end
 
 square1_boundary_model = rigid_body_boundary_model(square1)
@@ -234,15 +235,11 @@ contact_model = RigidContactModel(; normal_stiffness=2.0e5,
                                   contact_distance=2.0 * structure_particle_spacing)
 nothing # hide
 
-rigid_body_system_1_step3 = RigidBodySystem(square1;
-                                            boundary_model=square1_boundary_model,
-                                            contact_model=contact_model,
-                                            acceleration=(0.0, -gravity),
+rigid_body_system_1_step3 = RigidBodySystem(square1; boundary_model=square1_boundary_model,
+                                            contact_model, acceleration=(0.0, -gravity),
                                             particle_spacing=structure_particle_spacing)
-rigid_body_system_2_step3 = RigidBodySystem(square2;
-                                            boundary_model=square2_boundary_model,
-                                            contact_model=contact_model,
-                                            acceleration=(0.0, -gravity),
+rigid_body_system_2_step3 = RigidBodySystem(square2; boundary_model=square2_boundary_model,
+                                            contact_model, acceleration=(0.0, -gravity),
                                             particle_spacing=structure_particle_spacing)
 nothing # hide
 
@@ -301,15 +298,11 @@ circle1_boundary_model = rigid_body_boundary_model(circle1)
 circle2_boundary_model = rigid_body_boundary_model(circle2)
 nothing # hide
 
-rigid_body_system_1_step4 = RigidBodySystem(circle1;
-                                            boundary_model=circle1_boundary_model,
-                                            contact_model=contact_model,
-                                            acceleration=(0.0, -gravity),
+rigid_body_system_1_step4 = RigidBodySystem(circle1; boundary_model=circle1_boundary_model,
+                                            contact_model, acceleration=(0.0, -gravity),
                                             particle_spacing=structure_particle_spacing)
-rigid_body_system_2_step4 = RigidBodySystem(circle2;
-                                            boundary_model=circle2_boundary_model,
-                                            contact_model=contact_model,
-                                            acceleration=(0.0, -gravity),
+rigid_body_system_2_step4 = RigidBodySystem(circle2; boundary_model=circle2_boundary_model,
+                                            contact_model, acceleration=(0.0, -gravity),
                                             particle_spacing=structure_particle_spacing)
 nothing # hide
 
@@ -359,10 +352,8 @@ small_spheres = [SphereShape(structure_particle_spacing, small_sphere_radius,
 
 small_sphere_systems = [begin
                             sphere_boundary_model = rigid_body_boundary_model(sphere)
-                            RigidBodySystem(sphere;
-                                            boundary_model=sphere_boundary_model,
-                                            contact_model=contact_model,
-                                            acceleration=(0.0, -gravity),
+                            RigidBodySystem(sphere; boundary_model=sphere_boundary_model,
+                                            contact_model, acceleration=(0.0, -gravity),
                                             particle_spacing=structure_particle_spacing)
                         end
                         for sphere in small_spheres]
@@ -418,10 +409,8 @@ hexagon_shape = TrixiParticles.@set hexagon_shape.particle_spacing = structure_p
 # just like the other shapes in this tutorial.
 hexagon_boundary_model = rigid_body_boundary_model(hexagon_shape)
 
-hexagon_system = RigidBodySystem(hexagon_shape;
-                                 boundary_model=hexagon_boundary_model,
-                                 contact_model=contact_model,
-                                 acceleration=(0.0, -gravity),
+hexagon_system = RigidBodySystem(hexagon_shape; boundary_model=hexagon_boundary_model,
+                                 contact_model, acceleration=(0.0, -gravity),
                                  particle_spacing=structure_particle_spacing)
 #
 # You can then create a semidiscretization with this new system.