Implement OscillatingMotion (#915)

* Add `OscillatingMotion`

* Fix docs

* Add tests

* Reformat

* Fix docs

* Change API to use keyword arguments

* Implement suggestions

---------

Co-authored-by: Niklas Neher <[email protected]>
Co-authored-by: Sven Berger <[email protected]>

Author: 3 people

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")]
 ```
 
 

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/schemes/boundary/prescribed_motion.jl

@@ -10,7 +10,7 @@
                        are moving or not. Its boolean return value determines
                        if the neighborhood search will be updated.
 
-# Keyword Arguments
+# Keywords
 - `moving_particles`: Indices of moving particles. Default is each particle in the system.
 
 # Examples
@@ -98,3 +98,83 @@ function (prescribed_motion::Nothing)(system::AbstractSystem, t, semi)
 
     return system
 end
+
+"""
+    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 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.

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

src/schemes/fluid/weakly_compressible_sph/system.jl

@@ -26,7 +26,7 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method.
 - `smoothing_length`:   Smoothing length to be used for this system.
                         See [Smoothing Kernels](@ref smoothing_kernel).
 
-# Keyword Arguments
+# Keywords
 - `acceleration`:               Acceleration vector for the system. (default: zero vector)
 - `viscosity`:                  Viscosity model for this system (default: no viscosity).
                                 See [`ArtificialViscosityMonaghan`](@ref) or [`ViscosityAdami`](@ref).

src/schemes/structure/total_lagrangian_sph/system.jl

@@ -22,7 +22,7 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method.
 - `smoothing_length`:   Smoothing length to be used for this system.
                         See [Smoothing Kernels](@ref smoothing_kernel).
 
-# Keyword Arguments
+# Keywords
 - `n_clamped_particles`: Number of clamped particles which are fixed and not integrated
                     to clamp the structure. Note that the clamped particles must be the **last**
                     particles in the `InitialCondition`. See the info box below.

test/schemes/boundary/prescribed_motion.jl

@@ -0,0 +1,140 @@
+# Note that more tests for `PrescribedMotion` are in `test/systems/boundary_system.jl`.
+@testset verbose=true "PrescribedMotion" begin
+    @testset "OscillatingMotion2D" begin
+        @testset "Simple Rotation" begin
+            frequency = 1.0
+            translation_vector = SVector(0.0, 0.0)
+            rotation_angle = pi
+            rotation_center = SVector(0.0, 0.0)
+
+            motion = OscillatingMotion2D(; frequency, translation_vector, rotation_angle,
+                                         rotation_center)
+            movement_function = motion.movement_function
+
+            @test isapprox(movement_function([1.0, 0.0], 0.0), [1.0, 0.0], atol=eps())
+            # Initial position after one phase
+            @test isapprox(movement_function([1.0, 0.0], 1.0), [1.0, 0.0], atol=10eps())
+            # Initial position after one half phase
+            @test isapprox(movement_function([1.0, 0.0], 0.5), [1.0, 0.0], atol=10eps())
+            # Half rotation after one quarter phase (rotation angle is one half rotation)
+            @test isapprox(movement_function([1.0, 0.0], 0.25), [-1.0, 0.0], atol=10eps())
+            # Quarter rotation (0.5 * rotation_angle = 90 degrees)
+            @test isapprox(movement_function([1.0, 0.0], asin(0.5) / 2pi), [0.0, 1.0],
+                           atol=10eps())
+        end
+
+        @testset "Simple Translation" begin
+            frequency = 1.0
+            translation_vector = SVector(1.0, 0.0)
+            rotation_angle = 0.0
+            rotation_center = SVector(0.0, 0.0)
+
+            motion = OscillatingMotion2D(; frequency, translation_vector, rotation_angle,
+                                         rotation_center)
+            movement_function = motion.movement_function
+
+            @test isapprox(movement_function([0.0, 0.0], 0.0), [0.0, 0.0], atol=eps())
+            @test isapprox(movement_function([0.0, 0.0], 0.25), [1.0, 0.0], atol=10eps())
+            @test isapprox(movement_function([0.0, 0.0], 0.5), [0.0, 0.0], atol=10eps())
+            @test isapprox(movement_function([0.0, 0.0], 0.75), [-1.0, 0.0], atol=10eps())
+            @test isapprox(movement_function([0.0, 0.0], 1.0), [0.0, 0.0], atol=10eps())
+        end
+
+        @testset "Simple Translation tspan" begin
+            frequency = 1.0
+            translation_vector = SVector(1.0, 0.0)
+            rotation_angle = 0.0
+            rotation_center = SVector(0.0, 0.0)
+            tspan = (0.4, 1.4)
+
+            motion = OscillatingMotion2D(; frequency, translation_vector, rotation_angle,
+                                         rotation_center, tspan)
+            movement_function = motion.movement_function
+
+            @test isapprox(movement_function([0.0, 0.0], 0.4), [0.0, 0.0], atol=eps())
+            @test isapprox(movement_function([0.0, 0.0], 0.65), [1.0, 0.0], atol=10eps())
+            @test isapprox(movement_function([0.0, 0.0], 0.9), [0.0, 0.0], atol=10eps())
+            @test isapprox(movement_function([0.0, 0.0], 1.15), [-1.0, 0.0], atol=10eps())
+            @test isapprox(movement_function([0.0, 0.0], 1.4), [0.0, 0.0], atol=10eps())
+        end
+
+        @testset "Combined Rotation and Translation" begin
+            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)
+            movement_function = motion.movement_function
+
+            @test isapprox(movement_function([1.0, 0.0], 0.0), [1.0, 0.0], atol=eps())
+            @test isapprox(movement_function([1.0, 0.0], 0.25), [1.0, 1.0], atol=10eps())
+            @test isapprox(movement_function([1.0, 0.0], 0.5), [1.0, 0.0], atol=10eps())
+            @test isapprox(movement_function([1.0, 0.0], 0.75), [-1.0, -1.0], atol=10eps())
+            @test isapprox(movement_function([1.0, 0.0], 1.0), [1.0, 0.0], atol=10eps())
+        end
+
+        @testset "Phase Offset" begin
+            frequency = 1.0
+            translation_vector = SVector(1.0, 0.0)
+            rotation_angle = pi / 2
+            rotation_center = SVector(0.0, 0.0)
+            rotation_phase_offset = 0.25
+
+            motion = OscillatingMotion2D(; frequency, translation_vector, rotation_angle,
+                                         rotation_center, rotation_phase_offset)
+            movement_function = motion.movement_function
+
+            @test isapprox(movement_function([1.0, 0.0], 0.0), [0.0, -1.0], atol=eps())
+            @test isapprox(movement_function([1.0, 0.0], 0.25), [2.0, 0.0], atol=10eps())
+            @test isapprox(movement_function([1.0, 0.0], 0.5), [0.0, 1.0], atol=10eps())
+            @test isapprox(movement_function([1.0, 0.0], 0.75), [0.0, 0.0], atol=10eps())
+            @test isapprox(movement_function([1.0, 0.0], 1.0), [0.0, -1.0], atol=10eps())
+        end
+
+        @testset "Ramp-up" begin
+            frequency = 1.0
+            translation_vector = SVector(1.0, 0.0)
+            rotation_angle = 0.0
+            rotation_center = SVector(0.0, 0.0)
+            ramp_up_tspan = (0.0, 1.0)
+
+            motion = OscillatingMotion2D(; frequency, translation_vector, rotation_angle,
+                                         rotation_center, ramp_up_tspan)
+            movement_function = motion.movement_function
+
+            @test isapprox(movement_function([0.0, 0.0], 0.0), [0.0, 0.0], atol=eps())
+            # During ramp-up, amplitude is reduced
+            @test isapprox(movement_function([0.0, 0.0], 0.25), [0.15625, 0.0],
+                           atol=10eps())
+            @test isapprox(movement_function([0.0, 0.0], 0.5), [0.0, 0.0], atol=10eps())
+            # Less reduction at later times during ramp-up
+            @test isapprox(movement_function([0.0, 0.0], 0.75), [-0.84375, 0.0],
+                           atol=10eps())
+            @test isapprox(movement_function([0.0, 0.0], 1.0), [0.0, 0.0], atol=10eps())
+            # After ramp-up, full amplitude
+            @test isapprox(movement_function([0.0, 0.0], 1.25), [1.0, 0.0], atol=10eps())
+        end
+
+        @testset "tspan" begin
+            frequency = 1.0
+            translation_vector = SVector(1.0, 0.0)
+            rotation_angle = 0.0
+            rotation_center = SVector(0.0, 0.0)
+            tspan = (0.5, 1.5)
+
+            motion = OscillatingMotion2D(; frequency, translation_vector, rotation_angle,
+                                         rotation_center, tspan)
+            ismoving = motion.is_moving
+
+            @test !ismoving(-0.5)
+            @test !ismoving(0.0)
+            @test !ismoving(0.25)
+            @test ismoving(0.5)
+            @test ismoving(1.0)
+            @test ismoving(1.5)
+            @test !ismoving(1.75)
+        end
+    end
+end