Refactor SPH system initialization to use keyword arguments

NEWS.md

@@ -16,6 +16,53 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 - `DensityDiffusionAntuono` now only takes only one kwarg `delta` (#1142).
 - Return type of `vtk2trixi` changed to `NamedTuple` including an optional
   `:initial_condition` field if `create_initial_condition=true` is passed (#959).
+- Public system constructors now use keyword arguments for configuration values.
+  This affects `WeaklyCompressibleSPHSystem`, `EntropicallyDampedSPHSystem`,
+  `ImplicitIncompressibleSPHSystem`, `TotalLagrangianSPHSystem`, `DEMSystem`,
+  and `BoundaryDEMSystem`.
+
+#### Migration Guide for Keyword-Based System Constructors
+
+The initial condition remains the only positional argument for the affected system
+constructors. All method, material, and configuration arguments that were previously
+passed positionally now need to be passed as keywords.
+Existing optional keyword arguments remain available. For example, the
+`density_calculator` keyword of `EntropicallyDampedSPHSystem` still defaults to
+`SummationDensity()`.
+
+Here are examples, where `ic` is the initial condition, `kernel` is the
+smoothing kernel, and `h` is the smoothing length:
+
+```julia
+WeaklyCompressibleSPHSystem(ic, density_calculator, state_equation, kernel, h;
+                            kwargs...)
+# becomes
+WeaklyCompressibleSPHSystem(ic; smoothing_kernel=kernel, smoothing_length=h,
+                            density_calculator, state_equation, kwargs...)
+
+EntropicallyDampedSPHSystem(ic, kernel, h, sound_speed; kwargs...)
+# becomes
+EntropicallyDampedSPHSystem(ic; smoothing_kernel=kernel, smoothing_length=h,
+                            sound_speed, kwargs...)
+
+ImplicitIncompressibleSPHSystem(ic, kernel, h, reference_density; kwargs...)
+# becomes
+ImplicitIncompressibleSPHSystem(ic; smoothing_kernel=kernel, smoothing_length=h,
+                                reference_density, kwargs...)
+
+TotalLagrangianSPHSystem(ic, kernel, h, young_modulus, poisson_ratio; kwargs...)
+# becomes
+TotalLagrangianSPHSystem(ic; smoothing_kernel=kernel, smoothing_length=h,
+                         young_modulus, poisson_ratio, kwargs...)
+
+DEMSystem(ic, contact_model; kwargs...)
+# becomes
+DEMSystem(ic; contact_model, kwargs...)
+
+BoundaryDEMSystem(ic, normal_stiffness)
+# becomes
+BoundaryDEMSystem(ic; normal_stiffness)
+```
 
 ### Performance
 

README.md

@@ -134,7 +134,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

@@ -118,13 +118,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,
@@ -147,9 +147,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
@@ -164,8 +162,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)
@@ -177,11 +174,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)
@@ -193,8 +187,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)
@@ -215,12 +208,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)
@@ -232,8 +223,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)
@@ -255,11 +245,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.
@@ -274,10 +261,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,
@@ -294,8 +279,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)
 
@@ -324,24 +309,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)
 
@@ -352,8 +330,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)