Packing gpu

NEWS.md

@@ -4,14 +4,27 @@ 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.4
+
+### API Changes
+
+- Combined transport velocity formulation (TVF) and particle shifting technique (PST) into
+  one unified framework.
+  The keyword argument `transport_velocity` now changed to `shifting_technique`.
+  The `ParticleShiftingCallback` has been removed. To use PST, use the `UpdateCallback`
+  instead, and pass `shifting_technique=ParticleShiftingTechnique()` to the system.
+
+- Renamed the keyword argument `tlsph` to `place_on_shell` for `ParticlePackingSystem`,
+  `sample_boundary`, `extrude_geometry`, `RectangularShape`, and `SphereShape`.
+
 ## Version 0.3.1
 
 ### Features
 
-- **Simplified SGS Viscosity Models**: Added ViscosityMorrisSGS and ViscosityAdamiSGS, 
+- **Simplified SGS Viscosity Models**: Added ViscosityMorrisSGS and ViscosityAdamiSGS,
   which implement a simplified Smagorinsky-type sub-grid-scale viscosity. (#753)
 
-- **Multithreaded Integration Array**: Introduced a new array type for CPU backends 
+- **Multithreaded Integration Array**: Introduced a new array type for CPU backends
   that enables multithreaded broadcasting, delivering speed-ups of up to 5× on systems
   with many threads when combined with thread pinning. (#722)
 
@@ -21,17 +34,17 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 - **DXF file format support**: Import complex geometries using the DXF file format. (#821)
 
 - **Improved Plane interpolation**: Massively improved interpolation performance for planes (#763).
-  
+
 ### GPU
 
  - Make PST GPU-compatible (#813).
-   
+
  - Make open boundaries GPU-compatible (#773).
-  
+
  - Make interpolation GPU-compatible (#812).
 
 ### Important Bugfixes
- 
+
  - Fix validation setups (#801).
 
  - Calculate interpolated density instead of computed density when using interpolation (#808).

Project.toml

@@ -5,6 +5,7 @@ version = "0.3.2-dev"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+ArraysOfArrays = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -41,6 +42,7 @@ TrixiParticlesOrdinaryDiffEqExt = ["OrdinaryDiffEq", "OrdinaryDiffEqCore"]
 
 [compat]
 Adapt = "4"
+ArraysOfArrays = "0.6.5"
 CSV = "0.10"
 DataFrames = "1.6"
 DelimitedFiles = "1"

docs/src/systems/entropically_damped_sph.md

@@ -45,64 +45,3 @@ is a good choice for a wide range of Reynolds numbers (0.0125 to 10000).
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "fluid", "entropically_damped_sph", "system.jl")]
 ```
-
-## [Transport Velocity Formulation (TVF)](@id transport_velocity_formulation)
-Standard SPH suffers from problems like tensile instability or the creation of void regions in the flow.
-To address these problems, [Adami (2013)](@cite Adami2013) modified the advection velocity and added an extra term to the momentum equation.
-The authors introduced the so-called Transport Velocity Formulation (TVF) for WCSPH. [Ramachandran (2019)](@cite Ramachandran2019) applied the TVF
-also for the [EDAC](@ref edac) scheme.
-
-The transport velocity ``\tilde{v}_a`` of particle ``a`` is used to evolve the position of the particle ``r_a`` from one time step to the next by
-
-```math
-\frac{\mathrm{d} r_a}{\mathrm{d}t} = \tilde{v}_a
-```
-
-and is obtained at every time-step ``\Delta t`` from
-
-```math
-\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right),
-```
-
-where ``\rho_a`` is the density of particle ``a`` and ``p_{\text{background}}`` is a constant background pressure field.
-The tilde in the second term of the right hand side indicates that the material derivative has an advection part.
-
-The discretized form of the last term is
-
-```math
- -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx  -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab},
-```
-
-where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively.
-Note that although in the continuous case ``\nabla p_{\text{background}} = 0``, the discretization is not 0th-order consistent for **non**-uniform particle distribution,
-which means that there is a non-vanishing contribution only when particles are disordered.
-That also means that ``p_{\text{background}}`` occurs as prefactor to correct the trajectory of a particle resulting in uniform pressure distributions.
-Suggested is a background pressure which is in the order of the reference pressure but can be chosen arbitrarily large when the time-step criterion is adjusted.
-
-The inviscid momentum equation with an additional convection term for a particle moving with ``\tilde{v}`` is
-
-```math
-\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p +  \nabla \cdot \bm{A},
-```
-
- where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)^T`` is a consequence of the modified
- advection velocity and can be interpreted as the convection of momentum with the relative velocity ``\tilde{v}-v``.
-
-The discretized form of the momentum equation for a particle ``a`` reads as
-
-```math
-\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right].
-```
-
-Here, ``\tilde{p}_{ab}`` is the density-weighted pressure
-
-```math
-\tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b},
-```
-
-with the density  ``\rho_a``,  ``\rho_b`` and the pressure  ``p_a``,  ``p_b`` of particles ``a`` and ``b`` respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors for particle ``a`` and ``b`` respectively and are given, e.g. for particle ``a``, as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``.
-
-```@autodocs
-Modules = [TrixiParticles]
-Pages = [joinpath("schemes", "fluid", "transport_velocity.jl")]
-```

docs/src/systems/weakly_compressible_sph.md

@@ -152,8 +152,74 @@ as explained in [Sun2018](@cite Sun2018) on page 29, right above Equation 9.
 The ``\delta``-SPH method (WCSPH with density diffusion) together with this formulation
 of PST is commonly referred to as ``\delta^+``-SPH.
 
-The Particle Shifting Technique can be applied in form
-of the [`ParticleShiftingCallback`](@ref).
+To apply particle shifting, use the keyword argument `shifting_technique` in the constructor
+of a system that supports it.
+
+
+## [Transport Velocity Formulation (TVF)](@id transport_velocity_formulation)
+
+An alternative formulation is the so-called Transport Velocity Formulation (TVF)
+by [Adami (2013)](@cite Adami2013).
+[Ramachandran (2019)](@cite Ramachandran2019) applied the TVF also for the [EDAC](@ref edac)
+scheme.
+
+The transport velocity ``\tilde{v}_a`` of particle ``a`` is used to evolve the position
+of the particle ``r_a`` from one time step to the next by
+```math
+\frac{\mathrm{d} r_a}{\mathrm{d}t} = \tilde{v}_a
+```
+and is obtained at every time step ``\Delta t`` from
+```math
+\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right),
+```
+where ``\rho_a`` is the density of particle ``a`` and ``p_{\text{background}}``
+is a constant background pressure field.
+The tilde in the second term of the right-hand side indicates that the material derivative
+has an advection part.
+
+The discretized form of the last term is
+```math
+ -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx  -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab},
+```
+where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively.
+Note that although in the continuous case ``\nabla p_{\text{background}} = 0``,
+the discretization is not 0th-order consistent for **non**-uniform particle distribution,
+which means that there is a non-vanishing contribution only when particles are disordered.
+That also means that ``p_{\text{background}}`` occurs as pre-factor to correct
+the trajectory of a particle resulting in uniform pressure distributions.
+Suggested is a background pressure which is in the order of the reference pressure,
+but it can be chosen arbitrarily large when the time-step criterion is adjusted.
+
+The inviscid momentum equation with an additional convection term for a particle
+moving with ``\tilde{v}`` is
+```math
+\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p +  \nabla \cdot \bm{A},
+```
+where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)^T`` is a consequence
+of the modified advection velocity and can be interpreted as the convection of momentum
+with the relative velocity ``\tilde{v}-v``.
+
+The discretized form of the momentum equation for a particle ``a`` reads as
+```math
+\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right].
+```
+Here, ``\tilde{p}_{ab}`` is the density-weighted pressure
+```math
+\tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b},
+```
+with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a``
+and ``b``, respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors
+for particle ``a`` and ``b``, respectively, and are given, e.g., for particle ``a``,
+as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``.
+
+To apply the TVF, use the keyword argument `shifting_technique` in the constructor
+of a system that supports it.
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "fluid", "shifting_techniques.jl")]
+```
+
 
 ## [Tensile Instability Control](@id tic)
 

examples/dem/collapsing_sand_pile_3d.jl

@@ -55,7 +55,8 @@ min_coords_floor = (min_boundary[1] - boundary_thickness,
 floor_particles = RectangularShape(particle_spacing,
                                    (n_particles_floor_x, n_particles_floor_y,
                                     n_particles_floor_z),
-                                   min_coords_floor; density=boundary_density, tlsph=true)
+                                   min_coords_floor; density=boundary_density,
+                                   place_on_shell=true)
 boundary_particles = floor_particles
 
 # ==========================================================================================

examples/fluid/lid_driven_cavity_2d.jl

@@ -65,15 +65,15 @@ if wcsph
                                                state_equation, smoothing_kernel,
                                                pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure,
                                                smoothing_length, viscosity=viscosity,
-                                               transport_velocity=TransportVelocityAdami(pressure))
+                                               shifting_technique=TransportVelocityAdami(pressure))
 else
     state_equation = nothing
     density_calculator = ContinuityDensity()
     fluid_system = EntropicallyDampedSPHSystem(cavity.fluid, smoothing_kernel,
                                                smoothing_length,
                                                density_calculator=density_calculator,
                                                sound_speed, viscosity=viscosity,
-                                               transport_velocity=TransportVelocityAdami(pressure))
+                                               shifting_technique=TransportVelocityAdami(pressure))
 end
 
 # ==========================================================================================

examples/fluid/periodic_array_of_cylinders_2d.jl

@@ -60,7 +60,7 @@ smoothing_length = 1.2 * particle_spacing
 smoothing_kernel = SchoenbergQuarticSplineKernel{2}()
 fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length,
                                            sound_speed, viscosity=ViscosityAdami(; nu),
-                                           transport_velocity=TransportVelocityAdami(pressure),
+                                           shifting_technique=TransportVelocityAdami(pressure),
                                            acceleration=(acceleration_x, 0.0))
 
 # ==========================================================================================

examples/fluid/periodic_channel_2d.jl

@@ -28,7 +28,7 @@ initial_fluid_size = tank_size
 initial_velocity = (1.0, 0.0)
 
 fluid_density = 1000.0
-sound_speed = initial_velocity[1]
+sound_speed = 10 * initial_velocity[1]
 state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
                                    exponent=7)
 
@@ -48,6 +48,7 @@ viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
 fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                            state_equation, smoothing_kernel,
                                            smoothing_length, viscosity=viscosity,
+                                           shifting_technique=nothing,
                                            pressure_acceleration=nothing)
 
 # ==========================================================================================

examples/fluid/pipe_flow_2d.jl

@@ -74,6 +74,7 @@ viscosity = ViscosityAdami(nu=kinematic_viscosity)
 fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothing_length,
                                            sound_speed, viscosity=viscosity,
                                            density_calculator=fluid_density_calculator,
+                                           shifting_technique=ParticleShiftingTechnique(),
                                            buffer_size=n_buffer_particles)
 
 # Alternatively the WCSPH scheme can be used
@@ -86,6 +87,7 @@ if wcsph
     fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator,
                                                state_equation, smoothing_kernel,
                                                smoothing_length, viscosity=viscosity,
+                                               shifting_technique=ParticleShiftingTechnique(),
                                                buffer_size=n_buffer_particles)
 end
 
@@ -159,12 +161,10 @@ ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100)
 saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
-particle_shifting = ParticleShiftingCallback()
 
 extra_callback = nothing
 
-callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(),
-                        particle_shifting, extra_callback)
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), extra_callback)
 
 sol = solve(ode, RDPK3SpFSAL35(),
             abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)

examples/fluid/taylor_green_vortex_2d.jl

@@ -86,13 +86,13 @@ if wcsph
                                                pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure,
                                                smoothing_length,
                                                viscosity=ViscosityAdami(; nu),
-                                               transport_velocity=TransportVelocityAdami(background_pressure))
+                                               shifting_technique=TransportVelocityAdami(background_pressure))
 else
     density_calculator = SummationDensity()
     fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length,
                                                sound_speed,
                                                density_calculator=density_calculator,
-                                               transport_velocity=TransportVelocityAdami(background_pressure),
+                                               shifting_technique=TransportVelocityAdami(background_pressure),
                                                viscosity=ViscosityAdami(; nu))
 end
 

examples/fsi/dam_break_gate_2d.jl

@@ -83,18 +83,18 @@ solid_particle_spacing = thickness / (n_particles_x - 1)
 n_particles_y = round(Int, length_beam / solid_particle_spacing) + 1
 
 # The bottom layer is sampled separately below. Note that the `RectangularShape` puts the
-# first particle half a particle spacing away from the boundary, which is correct for fluids,
-# but not for solids. We therefore need to pass `tlsph=true`.
+# first particle half a particle spacing away from the shell of the shape, which is
+# correct for fluids, but not for solids. We therefore need to pass `place_on_shell=true`.
 #
 # The right end of the plate is 0.2 from the right end of the tank.
 plate_position = 0.6 - n_particles_x * solid_particle_spacing
 plate = RectangularShape(solid_particle_spacing,
                          (n_particles_x, n_particles_y - 1),
                          (plate_position, solid_particle_spacing),
-                         density=solid_density, tlsph=true)
+                         density=solid_density, place_on_shell=true)
 fixed_particles = RectangularShape(solid_particle_spacing,
                                    (n_particles_x, 1), (plate_position, 0.0),
-                                   density=solid_density, tlsph=true)
+                                   density=solid_density, place_on_shell=true)
 
 solid = union(plate, fixed_particles)
 

examples/fsi/dam_break_plate_2d.jl

@@ -57,15 +57,15 @@ solid_particle_spacing = thickness / (n_particles_x - 1)
 n_particles_y = round(Int, length_beam / solid_particle_spacing) + 1
 
 # The bottom layer is sampled separately below. Note that the `RectangularShape` puts the
-# first particle half a particle spacing away from the boundary, which is correct for fluids,
-# but not for solids. We therefore need to pass `tlsph=true`.
+# first particle half a particle spacing away from the shell of the shape, which is
+# correct for fluids, but not for solids. We therefore need to pass `place_on_shell=true`.
 plate = RectangularShape(solid_particle_spacing,
                          (n_particles_x, n_particles_y - 1),
                          (2initial_fluid_size[1], solid_particle_spacing),
-                         density=solid_density, tlsph=true)
+                         density=solid_density, place_on_shell=true)
 fixed_particles = RectangularShape(solid_particle_spacing,
                                    (n_particles_x, 1), (2initial_fluid_size[1], 0.0),
-                                   density=solid_density, tlsph=true)
+                                   density=solid_density, place_on_shell=true)
 
 solid = union(plate, fixed_particles)
 

examples/preprocessing/complex_shape_2d.jl

@@ -14,30 +14,28 @@
 using TrixiParticles
 using Plots
 
-particle_spacing = 0.05
+particle_spacing = 0.1
 
 filename = "inverted_open_curve"
-file = joinpath("examples", "preprocessing", "data", filename * ".asc")
+file = joinpath(examples_dir(), "preprocessing", "data", filename * ".asc")
 
-geometry = load_geometry(file)
+geometry = load_geometry(file; parallelization_backend=PolyesterBackend(),
+                         element_type=typeof(particle_spacing))
 
 trixi2vtk(geometry)
 
-point_in_geometry_algorithm = WindingNumberJacobson(; geometry,
+point_in_geometry_algorithm = WindingNumberJacobson(geometry;
+                                                    winding=HierarchicalWinding(geometry),
                                                     winding_number_factor=0.4,
-                                                    hierarchical_winding=true)
+                                                    store_winding_number=true)
 
 # Returns `InitialCondition`
 shape_sampled = ComplexShape(geometry; particle_spacing, density=1.0,
-                             store_winding_number=true,
                              point_in_geometry_algorithm)
 
-trixi2vtk(shape_sampled.initial_condition)
-
-coordinates = stack(shape_sampled.grid)
-# trixi2vtk(shape_sampled.signed_distance_field)
-# trixi2vtk(coordinates, w=shape_sampled.winding_numbers)
+trixi2vtk(shape_sampled)
 
 # Plot the winding number field
+coordinates = stack(point_in_geometry_algorithm.cache.grid)
 plot(InitialCondition(; coordinates, density=1.0, particle_spacing),
-     zcolor=shape_sampled.winding_numbers)
+     zcolor=point_in_geometry_algorithm.cache.winding_numbers)