Add ThrustCalculator

Author: efaulhaber

src/TrixiParticles.jl

@@ -74,6 +74,7 @@ export BoundaryZone, InFlow, OutFlow, BidirectionalFlow
 export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
        PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback,
        SplitIntegrationCallback, MechanicalWorkCalculator, calculated_mechanical_work,
+       ThrustCalculator, calculated_thrust,
        SortingCallback
 export ContinuityDensity, SummationDensity
 export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique,

src/general/mechanical_work_calculator.jl

@@ -212,3 +212,127 @@ function compute_structure_acceleration!(dv, system, eachparticle,
 
     return v
 end
+
+"""
+    ThrustCalculator(system::TotalLagrangianSPHSystem, semi;
+                     direction,
+                     eachparticle=eachparticle(system))
+
+Functor that computes the instantaneous hydrodynamic force exerted by the fluid on a
+[`TotalLagrangianSPHSystem`](@ref), projected onto `direction`.
+It can be passed as a custom quantity to [`PostprocessCallback`](@ref).
+
+For a fixed fluid-interacting structure in a channel flow, choose `direction` as the
+direction of useful force and multiply the recorded thrust by the corresponding reference
+speed to obtain useful mechanical power.
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
+
+# Arguments
+- `system`: The [`TotalLagrangianSPHSystem`](@ref) whose hydrodynamic force should be
+            monitored.
+- `semi`: The [`Semidiscretization`](@ref) that contains `system`.
+
+# Keywords
+- `direction`: Direction onto which the hydrodynamic force is projected. The vector is
+               normalized internally.
+- `eachparticle=eachparticle(system)`: Iterator selecting which particles contribute.
+
+# Examples
+```jldoctest; output = false, setup = :(system = TotalLagrangianSPHSystem(RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); smoothing_kernel=WendlandC2Kernel{2}(), smoothing_length=1.0, young_modulus=1.0, poisson_ratio=1.0); semi = (; systems=(system,), parallelization_backend=SerialBackend()))
+# Create a thrust calculator in x-direction
+thrust_calculator = ThrustCalculator(system, semi; direction=SVector(1.0, 0.0))
+
+# After postprocessing, retrieve the latest thrust value
+thrust = calculated_thrust(thrust_calculator)
+
+# output
+0.0
+```
+"""
+mutable struct ThrustCalculator{ELTYPE, DV, EP, D}
+    thrust       :: ELTYPE
+    system_index :: Int
+    dv           :: DV
+    eachparticle :: EP
+    direction    :: D
+end
+
+# This should dispatch on `TotalLagrangianSPHSystem`, but this name is not yet defined
+# due to the include order.
+function ThrustCalculator(system::AbstractStructureSystem, semi;
+                          direction, eachparticle=eachparticle(system))
+    ELTYPE = eltype(system)
+    system_index = system_indices(system, semi)
+
+    # Check vector length before converting to `SVector` to avoid extremely long
+    # compile times when accidentally passing large vectors.
+    if length(direction) != ndims(system)
+        throw(ArgumentError("length of `direction` must match the number of dimensions"))
+    end
+    direction_ = SVector(Tuple(direction))
+    if iszero(direction_)
+        throw(ArgumentError("`direction` must not be zero"))
+    end
+    direction_ = normalize(direction_)
+
+    # Allocate buffer to write hydrodynamic accelerations for all particles.
+    # `PostprocessCallback` transfers data to the CPU before calling custom quantities,
+    # so this can be a regular `Array` even when the simulation is running on a GPU.
+    dv = Array{ELTYPE, 2}(undef, (v_nvariables(system), nparticles(system)))
+
+    return ThrustCalculator(zero(ELTYPE), system_index, dv, eachparticle, direction_)
+end
+
+function reset!(calculator::ThrustCalculator)
+    calculator.thrust = zero(calculator.thrust)
+
+    return calculator
+end
+
+"""
+    calculated_thrust(calculator::ThrustCalculator)
+
+Get the latest projected hydrodynamic force from a [`ThrustCalculator`](@ref).
+"""
+function calculated_thrust(calculator::ThrustCalculator)
+    return calculator.thrust
+end
+
+function (calculator::ThrustCalculator)(system, dv_ode, du_ode, v_ode, u_ode, semi, t)
+    if system_indices(system, semi) != calculator.system_index
+        return nothing
+    end
+
+    thrust = update_thrust!(system, calculator.eachparticle, calculator.direction,
+                            calculator.dv, v_ode, u_ode, semi, t)
+
+    calculator.thrust = thrust
+
+    return calculator.thrust
+end
+
+function update_thrust!(system, eachparticle, direction, dv, v_ode, u_ode, semi, t)
+    # Note that the systems and NHS have already been updated by the
+    # `PostprocessCallback` before calling this function.
+    @trixi_timeit timer() "calculate thrust" begin
+        compute_structure_acceleration!(dv, system, eachparticle, true,
+                                        v_ode, u_ode, semi, t)
+
+        return projected_force(dv, system, eachparticle, direction)
+    end
+end
+
+function projected_force(dv, system, eachparticle, direction)
+    force = zero(eltype(system))
+
+    # Note that this is a reduction, so we cannot use `@threaded` here.
+    for particle in eachparticle
+        dv_particle = extract_svector(dv, system, particle)
+        force_particle = system.mass[particle] * dv_particle
+        force += dot(force_particle, direction)
+    end
+
+    return force
+end

test/callbacks/mechanical_work_calculator.jl

@@ -2,8 +2,12 @@
     # Mock system
     struct MockSystem <: TrixiParticles.AbstractStructureSystem{2}
         eltype::Type
+        mass::Vector
+
+        function MockSystem(ELTYPE)
+            new(ELTYPE, ELTYPE[1, 2, 3, 4])
+        end
     end
-    TrixiParticles.nparticles(::MockSystem) = 4
     Base.eltype(system::MockSystem) = system.eltype
 
     function TrixiParticles.create_neighborhood_search(neighborhood_search,
@@ -40,6 +44,100 @@
         @test eltype(calculator.t) == Float32
     end
 
+    @testset "ThrustCalculator constructor and force projection" begin
+        @test_throws UndefKeywordError ThrustCalculator(system64, semi64)
+
+        calculator = ThrustCalculator(system64, semi64; direction=SVector(1.0, 0.0))
+        @test calculator.system_index == 1
+        @test calculator.thrust == 0.0
+        @test calculator.dv isa Array{Float64, 2}
+        @test size(calculator.dv) == (2, 4)
+        @test calculator.eachparticle == eachparticle(system64)
+        @test calculator.direction == SVector(1.0, 0.0)
+        @test calculated_thrust(calculator) == 0.0
+
+        calculator = ThrustCalculator(system32, semi32; direction=(0.0, 2.0),
+                                      eachparticle=2:3)
+        @test eltype(calculator.thrust) == Float32
+        @test calculator.direction == SVector(0.0f0, 1.0f0)
+        @test calculator.eachparticle == 2:3
+
+        @test_throws ArgumentError ThrustCalculator(system64, semi64; direction=(0.0, 0.0))
+
+        dv = [2.0 -1.0 0.0 3.0
+              4.0 5.0 -2.0 1.0]
+        @test TrixiParticles.projected_force(dv, system64, eachparticle(system64),
+                                             SVector(1.0, 0.0)) == 12.0
+        @test TrixiParticles.projected_force(dv, system64, 2:3,
+                                             SVector(0.0, 1.0)) == 4.0
+
+        TrixiParticles.reset!(calculator)
+        @test calculated_thrust(calculator) == 0.0f0
+    end
+
+    @testset "ThrustCalculator FSI force" begin
+        particle_spacing = 1.0
+        smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+        smoothing_length = 1.0
+        fluid_density = 1000.0
+        structure_density = 2000.0
+        particle_volume = particle_spacing^2
+
+        state_equation = StateEquationCole(sound_speed=10.0,
+                                           reference_density=fluid_density,
+                                           exponent=1.0)
+
+        fluid_ic = InitialCondition(; coordinates=reshape([0.0, 0.0], 2, 1),
+                                    velocity=zeros(2, 1),
+                                    mass=[particle_volume * fluid_density],
+                                    density=[fluid_density], particle_spacing)
+
+        fluid_system = WeaklyCompressibleSPHSystem(fluid_ic; smoothing_kernel,
+                                                   smoothing_length,
+                                                   density_calculator=SummationDensity(),
+                                                   state_equation,
+                                                   reference_particle_spacing=particle_spacing)
+
+        structure_coordinates = reshape([1.5, 0.0], 2, 1)
+        structure_ic = InitialCondition(; coordinates=structure_coordinates,
+                                        velocity=zeros(2, 1),
+                                        mass=[particle_volume * structure_density],
+                                        density=[structure_density], particle_spacing)
+
+        boundary_model = BoundaryModelDummyParticles([fluid_density],
+                                                     [particle_volume * fluid_density],
+                                                     AdamiPressureExtrapolation(),
+                                                     smoothing_kernel, smoothing_length;
+                                                     state_equation,
+                                                     reference_particle_spacing=particle_spacing)
+
+        structure_system = TotalLagrangianSPHSystem(structure_ic; smoothing_kernel,
+                                                    smoothing_length,
+                                                    young_modulus=1.0e5,
+                                                    poisson_ratio=0.3,
+                                                    boundary_model)
+
+        semi_ = Semidiscretization(fluid_system, structure_system)
+        ode = semidiscretize(semi_, (0.0, 0.01))
+        semi = ode.p.semi
+
+        v_ode, u_ode = ode.u0.x
+        dv_ode = zero(v_ode)
+        TrixiParticles.kick!(dv_ode, v_ode, u_ode, ode.p, 0.0)
+
+        fluid = semi.systems[1]
+        structure = semi.systems[2]
+        dv_fluid = TrixiParticles.wrap_v(dv_ode, fluid, semi)
+
+        thrust = ThrustCalculator(structure, semi; direction=SVector(1.0, 0.0))
+        thrust(structure, dv_ode, nothing, v_ode, u_ode, semi, 0.0)
+
+        expected_force = -fluid.mass[1] * dv_fluid[1, 1]
+        @test !iszero(expected_force)
+        @test isapprox(calculated_thrust(thrust), expected_force;
+                       rtol=sqrt(eps()), atol=sqrt(eps()))
+    end
+
     @testset "update_mechanical_work_calculator!" begin
         # In the first test, we just move the 2x2 grid of particles up against gravity
         # and test that the accumulated work is just the potential energy difference.

test/examples/examples.jl

@@ -207,8 +207,9 @@
             mechanical_work1 = MechanicalWorkCalculator(tlsph_system_new, semi)
             mechanical_work2 = MechanicalWorkCalculator(tlsph_system_new, semi;
                                                         only_compute_force_on_fluid=true)
+            thrust = ThrustCalculator(tlsph_system_new, semi; direction=SVector(0.0, 1.0))
             postprocess_callback = PostprocessCallback(; interval=1, mechanical_work1,
-                                                       mechanical_work2,
+                                                       mechanical_work2, thrust,
                                                        write_file_interval=0)
 
             sol = @trixi_test_nowarn solve(ode, RDPK3SpFSAL35(), save_everystep=false,
@@ -233,6 +234,7 @@
             @test isapprox(calculated_mechanical_work(mechanical_work2),
                            expected_energy_fluid,
                            rtol=5e-4)
+            @test isfinite(calculated_thrust(thrust))
         end
 
         @trixi_testset "fsi/falling_water_column_2d.jl" begin