Combine PST and TVF into a unified framework

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).

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/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
 

src/TrixiParticles.jl

@@ -65,10 +65,9 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian
        BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem
 export BoundaryZone, InFlow, OutFlow, BidirectionalFlow
 export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
-       PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback,
-       ParticleShiftingCallback
+       PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback
 export ContinuityDensity, SummationDensity
-export PenaltyForceGanzenmueller, TransportVelocityAdami
+export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique
 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
        SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel,
        WendlandC6Kernel, SpikyKernel, Poly6Kernel

src/callbacks/callbacks.jl

@@ -32,4 +32,3 @@ include("post_process.jl")
 include("stepsize.jl")
 include("update.jl")
 include("steady_state_reached.jl")
-include("particle_shifting.jl")