Merge branch 'main' into oriented-bbox

Author: LasNikas

NEWS.md

@@ -4,6 +4,13 @@ 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.2
+
+### Features
+
+- Added `OscillatingMotion2D` to create an oscillating `PrescribedMotion` combining
+  translation and rotation (#915).
+
 ## Version 0.4.1
 
 ### Features

docs/src/systems/boundary.md

@@ -8,8 +8,9 @@
     BoundaryDEMSystem
 ```
 
-```@docs
-    PrescribedMotion
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "boundary", "prescribed_motion.jl")]
 ```
 
 

examples/structure/oscillating_beam_2d.jl

@@ -63,7 +63,8 @@ structure_system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothi
                                             material.E, material.nu,
                                             n_clamped_particles=nparticles(clamped_particles),
                                             acceleration=(0.0, -gravity),
-                                            penalty_force=nothing, viscosity=nothing)
+                                            penalty_force=nothing, viscosity=nothing,
+                                            clamped_particles_motion=nothing)
 
 # ==========================================================================================
 # ==== Simulation

src/TrixiParticles.jl

@@ -89,7 +89,7 @@ export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureEx
        BernoulliPressureExtrapolation
 export FirstOrderMirroring, ZerothOrderMirroring, SimpleMirroring
 export HertzContactModel, LinearContactModel
-export PrescribedMotion
+export PrescribedMotion, OscillatingMotion2D
 export examples_dir, validation_dir
 export trixi2vtk, vtk2trixi
 export RectangularTank, RectangularShape, SphereShape, ComplexShape

src/general/semidiscretization.jl

@@ -402,10 +402,14 @@ end
 function initialize_neighborhood_searches!(semi)
     foreach_system(semi) do system
         foreach_system(semi) do neighbor
+            # TODO Initialize after adapting to the GPU.
+            # Currently, this cannot use `semi.parallelization_backend`
+            # because data is still on the CPU.
             PointNeighbors.initialize!(get_neighborhood_search(system, neighbor, semi),
                                        initial_coordinates(system),
                                        initial_coordinates(neighbor),
-                                       eachindex_y=each_active_particle(neighbor))
+                                       eachindex_y=each_active_particle(neighbor),
+                                       parallelization_backend=PolyesterBackend())
         end
     end
 

src/io/write_vtk.jl

@@ -381,6 +381,9 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem)
                            initial_coords(system, particle)
                            for particle in eachparticle(system)]
 
+    vtk["is_clamped"] = vcat(fill(0, system.n_integrated_particles),
+                             fill(1, nparticles(system) - system.n_integrated_particles))
+
     vtk["lame_lambda"] = system.lame_lambda
     vtk["lame_mu"] = system.lame_mu
     vtk["young_modulus"] = system.young_modulus

src/schemes/boundary/prescribed_motion.jl

@@ -11,7 +11,9 @@
                        if the neighborhood search will be updated.
 
 # Keyword Arguments
-- `moving_particles`: Indices of moving particles. Default is each particle in the system.
+- `moving_particles`: Indices of moving particles. Default is each particle in the system
+                      for the [`WallBoundarySystem`](@ref) and all clamped particles
+                      for the [`TotalLagrangianSPHSystem`](@ref).
 
 # Examples
 ```jldoctest; output = false
@@ -48,7 +50,9 @@ function PrescribedMotion(movement_function, is_moving; moving_particles=nothing
     return PrescribedMotion(movement_function, is_moving, moving_particles)
 end
 
-function initialize!(prescribed_motion::PrescribedMotion, initial_condition)
+function initialize_prescribed_motion!(prescribed_motion::PrescribedMotion,
+                                       initial_condition,
+                                       n_clamped_particles=nparticles(initial_condition))
     # Test `movement_function` return type
     pos = extract_svector(initial_condition.coordinates,
                           Val(size(initial_condition.coordinates, 1)), 1)
@@ -57,23 +61,34 @@ function initialize!(prescribed_motion::PrescribedMotion, initial_condition)
               "Returning regular `Vector`s causes allocations and significant performance overhead."
     end
 
-    # Empty `moving_particles` means all particles are moving
+    # Empty `moving_particles` means all clamped particles are moving.
+    # For boundaries, all particles are considered clamped.
     if isempty(prescribed_motion.moving_particles)
         # Default is an empty vector, since the number of particles is not known when
         # instantiating `PrescribedMotion`.
-        resize!(prescribed_motion.moving_particles, nparticles(initial_condition))
-        prescribed_motion.moving_particles .= collect(1:nparticles(initial_condition))
+        resize!(prescribed_motion.moving_particles, n_clamped_particles)
+
+        # Clamped particles for TLSPH are the last `n_clamped_particles` in the system
+        first_particle = nparticles(initial_condition) - n_clamped_particles + 1
+        clamped_particles = first_particle:nparticles(initial_condition)
+        prescribed_motion.moving_particles .= collect(clamped_particles)
     end
+
+    return prescribed_motion
+end
+
+function initialize_prescribed_motion!(::Nothing, initial_condition,
+                                       n_clamped_particles=nparticles(initial_condition))
+    return nothing
 end
 
-function (prescribed_motion::PrescribedMotion)(system, t, semi)
-    (; coordinates, cache) = system
+function (prescribed_motion::PrescribedMotion)(coordinates, velocity, acceleration,
+                                               ismoving, system, semi, t)
     (; movement_function, is_moving, moving_particles) = prescribed_motion
-    (; acceleration, velocity) = cache
 
-    system.ismoving[] = is_moving(t)
+    ismoving[] = is_moving(t)
 
-    is_moving(t) || return system
+    is_moving(t) || return nothing
 
     @threaded semi for particle in moving_particles
         pos_original = initial_coords(system, particle)
@@ -90,11 +105,85 @@ function (prescribed_motion::PrescribedMotion)(system, t, semi)
         end
     end
 
-    return system
+    return nothing
 end
 
-function (prescribed_motion::Nothing)(system::AbstractSystem, t, semi)
-    system.ismoving[] = false
+"""
+    OscillatingMotion2D(; frequency, translation_vector, rotation_angle, rotation_center,
+                        rotation_phase_offset=0, tspan=(-Inf, Inf),
+                        ramp_up_tspan=(0.0, 0.0), moving_particles=nothing)
+
+Create a [`PrescribedMotion`](@ref) for a 2D oscillating motion of particles.
+The motion is a combination of a translation and a rotation around a center point
+with the same frequency. Note that a phase offset can be specified for a rotation
+that is out of sync with the translation.
+
+Create a [`PrescribedMotion`](@ref) for 2D oscillatory particle motion
+that combines a translation and a rotation about a specified center.
+Both use the same frequency; an optional phase offset allows the rotation
+to be out of phase with the translation.
+
+# Keywords
+- `frequency`:          Frequency of the oscillation in Hz.
+- `translation_vector`: Vector specifying the amplitude and direction of the translation.
+- `rotation_angle`:     Maximum rotation angle in radians.
+- `rotation_center`:    Center point of the rotation as an `SVector`.
+
+# Keywords (optional)
+- `rotation_phase_offset=0`: Phase offset of the rotation in number of periods (`1` = 1 period).
+- `tspan=(-Inf, Inf)`:  Time span in which the motion is active. Outside of this time span,
+                        particles remain at their last position.
+- `ramp_up_tspan`:      Time span in which the motion is smoothly ramped up from zero
+                        to full amplitude.
+                        Outside of this time span, the motion has full amplitude.
+                        Default is no ramp-up.
+- `moving_particles`:   Indices of moving particles. Default is each particle in the system
+                        for the [`WallBoundarySystem`](@ref) and all clamped particles
+                        for the [`TotalLagrangianSPHSystem`](@ref).
+
+# Examples
+```jldoctest; output = false, filter = r"PrescribedMotion{.*"
+# Oscillating motion with frequency 1 Hz, translation amplitude 1 in x-direction,
+# rotation angle 90 degrees (pi/2) around the origin.
+frequency = 1.0
+translation_vector = SVector(1.0, 0.0)
+rotation_angle = pi / 2
+rotation_center = SVector(0.0, 0.0)
+
+motion = OscillatingMotion2D(frequency, translation_vector, rotation_angle,
+                             rotation_center)
+
+# output
+PrescribedMotion{...}
+"""
+function OscillatingMotion2D(; frequency, translation_vector, rotation_angle,
+                             rotation_center, rotation_phase_offset=0, tspan=(-Inf, Inf),
+                             ramp_up_tspan=(0.0, 0.0), moving_particles=nothing)
+    @inline function movement_function(x, t)
+        if isfinite(tspan[1])
+            t = t - tspan[1]
+        end
+
+        sin_scaled = sin(frequency * 2pi * t)
+        translation = sin_scaled * translation_vector
+        x_centered = x .- rotation_center
+        sin_scaled_offset = sin(2pi * (frequency * t - rotation_phase_offset))
+        angle = rotation_angle * sin_scaled_offset
+        rotated = SVector(x_centered[1] * cos(angle) - x_centered[2] * sin(angle),
+                          x_centered[1] * sin(angle) + x_centered[2] * cos(angle))
+
+        result = rotated .+ rotation_center .+ translation
+
+        if ramp_up_tspan[2] > ramp_up_tspan[1] && ramp_up_tspan[1] <= t <= ramp_up_tspan[2]
+            # Apply a smoothstep ramp-up during the ramp-up time span
+            t_rel = (t - ramp_up_tspan[1]) / (ramp_up_tspan[2] - ramp_up_tspan[1])
+            ramp_factor = 3 * t_rel^2 - 2 * t_rel^3
+            return result * ramp_factor + (1 - ramp_factor) * x
+        end
+        return result
+    end
+
+    is_moving(t) = tspan[1] <= t <= tspan[2]
 
-    return system
+    return PrescribedMotion(movement_function, is_moving; moving_particles)
 end

src/schemes/boundary/wall_boundary/system.jl

@@ -9,7 +9,7 @@ The interaction between fluid and boundary particles is specified by the boundar
 - `initial_condition`: Initial condition (see [`InitialCondition`](@ref))
 - `boundary_model`: Boundary model (see [Boundary Models](@ref boundary_models))
 
-# Keyword Arguments
+# Keywords
 - `prescribed_motion`: For moving boundaries, a [`PrescribedMotion`](@ref) can be passed.
 - `adhesion_coefficient`: Coefficient specifying the adhesion of a fluid to the surface.
    Note: currently it is assumed that all fluids have the same adhesion coefficient.
@@ -42,12 +42,11 @@ function WallBoundarySystem(initial_condition, model; prescribed_motion=nothing,
     coordinates = copy(initial_condition.coordinates)
 
     ismoving = Ref(!isnothing(prescribed_motion))
-    initialize!(prescribed_motion, initial_condition)
+    initialize_prescribed_motion!(prescribed_motion, initial_condition)
 
     cache = create_cache_boundary(prescribed_motion, initial_condition)
     cache = (cache..., color=Int(color_value))
 
-    # Because of dispatches boundary model needs to be first!
     return WallBoundarySystem(initial_condition, coordinates, model, prescribed_motion,
                               ismoving, adhesion_coefficient, cache)
 end
@@ -105,7 +104,7 @@ end
     return current_velocity(v, system, system.prescribed_motion, particle)
 end
 
-@inline function current_velocity(v, system, movement, particle)
+@inline function current_velocity(v, system, prescribed_motion, particle)
     (; cache, ismoving) = system
 
     if ismoving[]
@@ -115,7 +114,7 @@ end
     return zero(SVector{ndims(system), eltype(system)})
 end
 
-@inline function current_velocity(v, system, movement::Nothing, particle)
+@inline function current_velocity(v, system, prescribed_motion::Nothing, particle)
     return zero(SVector{ndims(system), eltype(system)})
 end
 
@@ -127,7 +126,8 @@ end
     return current_acceleration(system, system.prescribed_motion, particle)
 end
 
-@inline function current_acceleration(system, ::PrescribedMotion, particle)
+@inline function current_acceleration(system::WallBoundarySystem, ::PrescribedMotion,
+                                      particle)
     (; cache, ismoving) = system
 
     if ismoving[]
@@ -137,7 +137,7 @@ end
     return zero(SVector{ndims(system), eltype(system)})
 end
 
-@inline function current_acceleration(system, ::Nothing, particle)
+@inline function current_acceleration(system::WallBoundarySystem, ::Nothing, particle)
     return zero(SVector{ndims(system), eltype(system)})
 end
 
@@ -177,7 +177,21 @@ end
 function update_positions!(system::WallBoundarySystem, v, u, v_ode, u_ode, semi, t)
     (; prescribed_motion) = system
 
-    prescribed_motion(system, t, semi)
+    apply_prescribed_motion!(system, prescribed_motion, semi, t)
+end
+
+function apply_prescribed_motion!(system::WallBoundarySystem,
+                                  prescribed_motion::PrescribedMotion, semi, t)
+    (; ismoving, coordinates, cache) = system
+    (; acceleration, velocity) = cache
+
+    prescribed_motion(coordinates, velocity, acceleration, ismoving, system, semi, t)
+
+    return system
+end
+
+function apply_prescribed_motion!(system::WallBoundarySystem, ::Nothing, semi, t)
+    return system
 end
 
 function update_quantities!(system::WallBoundarySystem, v, u, v_ode, u_ode, semi, t)

src/schemes/fluid/entropically_damped_sph/system.jl

@@ -22,7 +22,7 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more
 - `smoothing_length`:   Smoothing length to be used for this system.
                         See [Smoothing Kernels](@ref smoothing_kernel).
 
-# Keyword Arguments
+# Keywords
 - `viscosity`:                  Viscosity model for this system (default: no viscosity).
                                 Recommended: [`ViscosityAdami`](@ref).
 - `acceleration`:               Acceleration vector for the system. (default: zero vector)

src/schemes/fluid/surface_normal_sph.jl

@@ -4,7 +4,7 @@
 
 Color field based computation of the interface normals.
 
-# Keyword Arguments
+# Keywords
 - `boundary_contact_threshold=0.1`: If this threshold is reached the fluid is assumed to be in contact with the boundary.
 - `interface_threshold=0.01`:       Threshold for normals to be removed as being invalid.
 - `ideal_density_threshold=0.0`:    Assume particles are inside if they are above this threshold, which is relative to the `ideal_neighbor_count`.