diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 302657a193..50c151c22f 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -70,7 +70,7 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian export BoundaryZone, InFlow, OutFlow, BidirectionalFlow export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback, - SplitIntegrationCallback + SplitIntegrationCallback, EnergyCalculatorCallback, calculated_energy export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique, ParticleShiftingTechniqueSun2017, ConsistentShiftingSun2019, diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index b1bd602269..1be909f86a 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -33,3 +33,4 @@ include("stepsize.jl") include("update.jl") include("steady_state_reached.jl") include("split_integration.jl") +include("energy_calculator.jl") diff --git a/src/callbacks/energy_calculator.jl b/src/callbacks/energy_calculator.jl new file mode 100644 index 0000000000..9f5e9c917b --- /dev/null +++ b/src/callbacks/energy_calculator.jl @@ -0,0 +1,216 @@ +""" + EnergyCalculatorCallback(system, semi; interval=1, + eachparticle=(n_integrated_particles(system) + 1):nparticles(system), + only_compute_force_on_fluid=false) + +Callback to integrate the work done by particles in a +[`TotalLagrangianSPHSystem`](@ref). + +With the default arguments it tracks the energy contribution of the clamped particles +that follow a [`PrescribedMotion`](@ref). By selecting different particles it can also be +used to measure the work done by the structure on the surrounding fluid. + +- **Prescribed/clamped motion energy** (default) -- monitor only the clamped particles by + leaving `eachparticle` at its default range + `(n_integrated_particles(system) + 1):nparticles(system)`. +- **Fluid load measurement** -- set `eachparticle=eachparticle(system)` together with + `only_compute_force_on_fluid=true` to accumulate the work that the entire structure + exerts on the surrounding fluid (useful for drag or lift estimates). + +Internally the callback integrates the instantaneous power, i.e. the dot product between +the force exerted by the particle and its prescribed velocity, using an explicit Euler +time integration scheme. + +The accumulated value can be retrieved via [`calculated_energy`](@ref). + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. + +# Arguments +- `system`: The [`TotalLagrangianSPHSystem`](@ref) whose particles should be monitored. +- `semi`: The [`Semidiscretization`](@ref) that contains `system`. + +# Keywords +- `interval=1`: Interval (in number of time steps) at which to compute the energy. + It is recommended to keep this at `1` (every time step) or small (≤ 5) + to limit time integration errors in the energy integral. +- `eachparticle=(n_integrated_particles(system) + 1):nparticles(system)`: Iterator + selecting which particles contribute. The default includes all clamped + particles in the system; pass `eachparticle(system)` to include every particle. +- `only_compute_force_on_fluid=false`: When `true`, only interactions with + fluid systems are accounted for. Combined with + `eachparticle=eachparticle(system)`, this accumulates the work that the + entire structure exerts on the fluid, which is useful for drag or lift + estimates. + +# Examples +```jldoctest; output = false, setup = :(system = RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); semi = (; systems=(system,), parallelization_backend=SerialBackend())) +# Create an energy calculator callback that is called every 2 time steps +energy_cb = EnergyCalculatorCallback(system, semi; interval=2) + +# After the simulation, retrieve the calculated energy +total_energy = calculated_energy(energy_cb) + +# output +0.0 +``` +""" +struct EnergyCalculatorCallback{T, DV, EP} + interval :: Int + t :: T # Time of last call + energy :: T + system_index :: Int + dv :: DV + eachparticle :: EP + only_compute_force_on_fluid :: Bool +end + +function EnergyCalculatorCallback(system, semi; interval=1, + eachparticle=(n_integrated_particles(system) + 1):nparticles(system), + only_compute_force_on_fluid=false) + ELTYPE = eltype(system) + system_index = system_indices(system, semi) + + # Allocate buffer to write accelerations for all particles (including clamped ones) + dv = allocate(semi.parallelization_backend, ELTYPE, (ndims(system), nparticles(system))) + + cb = EnergyCalculatorCallback(interval, Ref(zero(ELTYPE)), Ref(zero(ELTYPE)), + system_index, dv, eachparticle, + only_compute_force_on_fluid) + + # The first one is the `condition`, the second the `affect!` + return DiscreteCallback(cb, cb, save_positions=(false, false)) +end + +""" + calculated_energy(cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback}) + +Get the current calculated energy from an [`EnergyCalculatorCallback`](@ref). + +# Arguments +- `cb`: The `DiscreteCallback` returned by [`EnergyCalculatorCallback`](@ref). + +# Examples +```jldoctest; output = false, setup = :(system = RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); semi = (; systems=(system,), parallelization_backend=SerialBackend())) +# Create an energy calculator callback +energy_cb = EnergyCalculatorCallback(system, semi) + +# After the simulation, retrieve the calculated energy +total_energy = calculated_energy(energy_cb) + +# output +0.0 +``` +""" +function calculated_energy(cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback}) + return cb.affect!.energy[] +end + +# `condition` +function (callback::EnergyCalculatorCallback)(u, t, integrator) + (; interval) = callback + + return condition_integrator_interval(integrator, interval) +end + +# `affect!` +function (callback::EnergyCalculatorCallback)(integrator) + # Determine time step size as difference to last time this callback was called + t = integrator.t + dt = t - callback.t[] + + # Update time of last call + callback.t[] = t + + semi = integrator.p + v_ode, u_ode = integrator.u.x + energy = callback.energy + + system = semi.systems[callback.system_index] + update_energy_calculator!(energy, system, callback.eachparticle, + callback.only_compute_force_on_fluid, callback.dv, + v_ode, u_ode, semi, t, dt) + + # Tell OrdinaryDiffEq that `u` has not been modified + u_modified!(integrator, false) + + return integrator +end + +function update_energy_calculator!(energy, system, eachparticle, + only_compute_force_on_fluid, dv, + v_ode, u_ode, semi, t, dt) + @trixi_timeit timer() "calculate energy" begin + # Update quantities that are stored in the systems. These quantities (e.g. pressure) + # still have the values from the last stage of the previous step if not updated here. + @trixi_timeit timer() "update systems and nhs" begin + # Don't create sub-timers here to avoid cluttering the timer output + @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t) + end + + set_zero!(dv) + + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + foreach_system(semi) do neighbor_system + if only_compute_force_on_fluid && !(neighbor_system isa AbstractFluidSystem) + # Not a fluid system, ignore this system + return + end + + v_neighbor = wrap_v(v_ode, neighbor_system, semi) + u_neighbor = wrap_u(u_ode, neighbor_system, semi) + + interact!(dv, v, u, v_neighbor, u_neighbor, + system, neighbor_system, semi, + integrate_tlsph=true, # Required when using split integration + eachparticle=eachparticle) + end + + if !only_compute_force_on_fluid + @threaded semi for particle in eachparticle + add_acceleration!(dv, particle, system) + add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t) + end + end + + for particle in eachparticle + velocity = current_velocity(v, system, particle) + dv_particle = extract_svector(dv, system, particle) + + # The force on the clamped particle is mass times acceleration + F_particle = system.mass[particle] * dv_particle + + # To obtain energy, we need to integrate the instantaneous power. + # Instantaneous power is force applied BY the particle times its velocity. + # The force applied BY the particle is the negative of the force applied ON it. + energy[] += dot(-F_particle, velocity) * dt + end + end + + return energy +end + +function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback}) + @nospecialize cb # reduce precompilation time + + ELTYPE = eltype(cb.affect!.energy) + print(io, "EnergyCalculatorCallback{$ELTYPE}(interval=", cb.affect!.interval, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + update_cb = cb.affect! + ELTYPE = eltype(update_cb.energy) + setup = [ + "interval" => update_cb.interval + ] + summary_box(io, "EnergyCalculatorCallback{$ELTYPE}", setup) + end +end diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl index 3e7b5fce32..d09b090bd9 100644 --- a/src/callbacks/split_integration.jl +++ b/src/callbacks/split_integration.jl @@ -1,6 +1,6 @@ mutable struct SplitIntegrationCallback # Type-stability does not matter here, and it is impossible to implement because - # the integration will be created during the initialization. + # the integrator will be created during the initialization. integrator :: Any alg :: Any kwargs :: Any @@ -71,9 +71,10 @@ function initialize_split_integration!(cb, u, t, integrator) # Create split integrator with TLSPH systems only systems = filter(i -> i isa TotalLagrangianSPHSystem, semi.systems) - # These neighborhood searches are never used + # This neighborhood search is never used + neighborhood_search = TrivialNeighborhoodSearch{ndims(first(systems))}() semi_split = Semidiscretization(systems..., - neighborhood_search=TrivialNeighborhoodSearch{ndims(first(systems))}(), + neighborhood_search=neighborhood_search, parallelization_backend=semi.parallelization_backend) # Verify that a NHS implementation is used that does not require updates diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index 8eecb8380a..056191330e 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -31,7 +31,7 @@ function UpdateCallback(; interval::Integer=-1, dt=0.0) update_callback! = UpdateCallback(interval) if dt > 0 - # Add a `tstop` every `dt`, and save the final solution. + # Add a `tstop` every `dt` return PeriodicCallback(update_callback!, dt, initialize=(initial_update!), save_positions=(false, false)) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index ceba8b5bdf..4c3206528f 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -232,7 +232,7 @@ end @inline function system_indices(system, semi) # Note that this takes only about 5 ns, while mapping systems to indices with a `Dict` # is ~30x slower because `hash(::System)` is very slow. - index = findfirst(==(system), semi.systems) + index = findfirst(s -> s === system, semi.systems) if isnothing(index) throw(ArgumentError("system is not in the semidiscretization")) diff --git a/src/schemes/structure/total_lagrangian_sph/rhs.jl b/src/schemes/structure/total_lagrangian_sph/rhs.jl index 3edad47a55..498cc78e59 100644 --- a/src/schemes/structure/total_lagrangian_sph/rhs.jl +++ b/src/schemes/structure/total_lagrangian_sph/rhs.jl @@ -3,18 +3,21 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, particle_system::TotalLagrangianSPHSystem, neighbor_system::TotalLagrangianSPHSystem, semi; - integrate_tlsph=semi.integrate_tlsph[]) + integrate_tlsph=semi.integrate_tlsph[], + eachparticle=each_integrated_particle(particle_system)) # Different structures do not interact with each other (yet) particle_system === neighbor_system || return dv # Skip interaction if TLSPH systems are integrated separately integrate_tlsph || return dv - interact_structure_structure!(dv, v_particle_system, particle_system, semi) + interact_structure_structure!(dv, v_particle_system, particle_system, semi; + eachparticle) end # Function barrier without dispatch for unit testing -@inline function interact_structure_structure!(dv, v_system, system, semi) +@inline function interact_structure_structure!(dv, v_system, system, semi; + eachparticle=each_integrated_particle(system)) (; penalty_force) = system # Everything here is done in the initial coordinates @@ -23,10 +26,10 @@ end # Loop over all pairs of particles and neighbors within the kernel cutoff. # For structure-structure interaction, this has to happen in the initial coordinates. foreach_point_neighbor(system, system, system_coords, system_coords, semi; - points=each_integrated_particle(system)) do particle, neighbor, - initial_pos_diff, - initial_distance - # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details. + points=eachparticle) do particle, neighbor, + initial_pos_diff, initial_distance + # Only consider particles with a distance > 0. + # See `src/general/smoothing_kernels.jl` for more details. initial_distance^2 < eps(initial_smoothing_length(system)^2) && return rho_a = @inbounds system.material_density[particle] @@ -72,7 +75,8 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, particle_system::TotalLagrangianSPHSystem, neighbor_system::AbstractFluidSystem, semi; - integrate_tlsph=semi.integrate_tlsph[]) + integrate_tlsph=semi.integrate_tlsph[], + eachparticle=each_integrated_particle(particle_system)) # Skip interaction if TLSPH systems are integrated separately integrate_tlsph || return dv @@ -84,10 +88,8 @@ function interact!(dv, v_particle_system, u_particle_system, # Loop over all pairs of particles and neighbors within the kernel cutoff foreach_point_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords, semi; - points=each_integrated_particle(particle_system)) do particle, - neighbor, - pos_diff, - distance + points=eachparticle) do particle, neighbor, + pos_diff, distance # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details. distance^2 < eps(initial_smoothing_length(particle_system)^2) && return @@ -178,7 +180,8 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, particle_system::TotalLagrangianSPHSystem, neighbor_system::Union{WallBoundarySystem, OpenBoundarySystem}, semi; - integrate_tlsph=semi.integrate_tlsph[]) + integrate_tlsph=semi.integrate_tlsph[], + eachparticle=each_integrated_particle(particle_system)) # TODO continuity equation? return dv end diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index 3979093ea0..0900ac42b9 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -336,11 +336,15 @@ end # Reset deformation gradient set_zero!(deformation_grad) - # Loop over all pairs of particles and neighbors within the kernel cutoff. + # Loop over all pairs of particles and neighbors within the kernel cutoff initial_coords = initial_coordinates(system) foreach_point_neighbor(system, system, initial_coords, initial_coords, - semi) do particle, neighbor, initial_pos_diff, initial_distance - # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details. + semi, + parallelization_backend=SerialBackend()) do particle, neighbor, + initial_pos_diff, + initial_distance + # Only consider particles with a distance > 0. + # See `src/general/smoothing_kernels.jl` for more details. initial_distance^2 < eps(initial_smoothing_length(system)^2) && return volume = mass[neighbor] / material_density[neighbor] diff --git a/test/callbacks/callbacks.jl b/test/callbacks/callbacks.jl index ae6166cde8..9c1ea952e7 100644 --- a/test/callbacks/callbacks.jl +++ b/test/callbacks/callbacks.jl @@ -5,4 +5,5 @@ include("update.jl") include("solution_saving.jl") include("steady_state_reached.jl") + include("energy_calculator.jl") end diff --git a/test/callbacks/energy_calculator.jl b/test/callbacks/energy_calculator.jl new file mode 100644 index 0000000000..74be5dc1b5 --- /dev/null +++ b/test/callbacks/energy_calculator.jl @@ -0,0 +1,152 @@ +@testset verbose=true "EnergyCalculatorCallback" begin + # Mock system + struct MockSystem <: TrixiParticles.AbstractSystem{2} end + TrixiParticles.nparticles(::MockSystem) = 4 + + function TrixiParticles.create_neighborhood_search(neighborhood_search, + system::MockSystem, neighbor) + return nothing + end + + system = MockSystem() + semi = Semidiscretization(system) + + @testset "Constructor and Basic Properties" begin + # Test default constructor + callback = EnergyCalculatorCallback{Float64}(system, semi) + @test callback.affect!.system_index == 1 + @test callback.affect!.interval == 1 + @test callback.affect!.t[] == 0.0 + @test callback.affect!.energy[] == 0.0 + @test callback.affect!.dv isa Array{Float64, 2} + @test size(callback.affect!.dv) == (2, 4) + @test callback.affect!.eachparticle == 1:4 + @test calculated_energy(callback) == 0.0 + + # Test constructor with interval + callback = EnergyCalculatorCallback{Float64}(system, semi; interval=5) + @test callback.affect!.interval == 5 + + # Test with specific element type + callback = EnergyCalculatorCallback{Float32}(system, semi; interval=2) + @test eltype(callback.affect!.energy) == Float32 + @test eltype(callback.affect!.t) == Float32 + end + + @testset "show" begin + callback = EnergyCalculatorCallback{Float64}(system, semi; interval=10) + + # Test compact representation + show_compact = "EnergyCalculatorCallback{Float64}(interval=10)" + @test repr(callback) == show_compact + + # Test detailed representation - check against expected box format + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ EnergyCalculatorCallback{Float64} │ + │ ═════════════════════════════════ │ + │ interval: ……………………………………………………… 10 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", callback) == show_box + end + + @testset "update_energy_calculator!" begin + # In the first test, we just move the 2x2 grid of particles up against gravity + # and test that the energy calculated is just the potential energy difference. + # In the other tests, we clamp the top row of particles and offset them to create + # stress. We can then test how much energy is required to pull the particles further + # apart against the elastic forces. + n_clamped_particles = [4, 2, 2, 2] + E = [1e6, 1e1, 1e6, 1e8] + @testset "n_clamped_particles = $(n_clamped_particles[i]), E = $(E[i])" for i in + eachindex(E) + # Create a simple 2D system with 2x2 particles + coordinates = [0.0 1.0 0.0 1.0; 0.0 0.0 1.0 1.0] + mass = [1.0, 1.0, 1.0, 1.0] + density = [1000.0, 1000.0, 1000.0, 1000.0] + + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = sqrt(2) + young_modulus = E[i] + poisson_ratio = 0.4 + + initial_condition = InitialCondition(; coordinates, mass, density) + + function movement_function(initial_pos, t) + # Offset clamped particles to create stress + return initial_pos + SVector(0.0, t + 0.5) + end + is_moving(t) = true + prescribed_motion = PrescribedMotion(movement_function, is_moving) + + # Create TLSPH system with energy calculator support + system = TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, + smoothing_length, young_modulus, + poisson_ratio, + clamped_particles_motion=prescribed_motion, + n_clamped_particles=n_clamped_particles[i], + acceleration=(0.0, -2.0), + use_with_energy_calculator_callback=true) + + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 1.0)) + + # Create dummy ODE state vectors + v = zeros(2, TrixiParticles.n_integrated_particles(system)) + u = coordinates[:, 1:TrixiParticles.n_integrated_particles(system)] + v_ode = vec(v) + u_ode = vec(u) + + # Update system + TrixiParticles.initialize!(system, semi) + TrixiParticles.update_positions!(system, v, u, v_ode, u_ode, semi, 0.0) + TrixiParticles.update_quantities!(system, v, u, v_ode, u_ode, semi, 0.0) + + # Set up test parameters + energy1 = Ref(0.0) + dt1 = 0.1 + + # Test that energy is integrated, i.e., values of instantaneous power + # are accumulated over time. This initial energy should just be an offset. + # Also, half the step size means half the energy increase. + energy2 = Ref(1.0) + dt2 = 0.05 + + # Test `only_compute_force_on_fluid` + energy3 = Ref(0.0) + dt3 = 0.1 + + eachparticle = (TrixiParticles.n_integrated_particles(system) + 1):nparticles(system) + dv = zeros(2, nparticles(system)) + + TrixiParticles.update_energy_calculator!(energy1, system, eachparticle, false, + dv, v_ode, u_ode, semi, 0.0, dt1) + TrixiParticles.update_energy_calculator!(energy2, system, eachparticle, false, + dv, v_ode, u_ode, semi, 0.0, dt2) + TrixiParticles.update_energy_calculator!(energy3, system, eachparticle, true, + dv, v_ode, u_ode, semi, 0.0, dt3) + + if i == 1 + @test isapprox(energy1[], 0.8) + @test isapprox(energy2[], 1.0 + 0.4) + elseif i == 2 + # For very soft material, we can just pull up the top row of particles + # and the energy required is almost just the potential energy difference. + @test isapprox(energy1[], 0.4080357142857143) + @test isapprox(energy2[], 1.0 + 0.5 * 0.4080357142857143) + elseif i == 3 + # For a stiffer material, the stress from the offset creates larger forces + # pulling the clamped particles back down, so we need a lot of energy + # to pull the material apart. + @test isapprox(energy1[], 803.9714285714281) + @test isapprox(energy2[], 1.0 + 0.5 * 803.9714285714281) + elseif i == 4 + # For a very stiff material, the energy is even larger. + @test isapprox(energy1[], 80357.5428571428) + @test isapprox(energy2[], 1.0 + 0.5 * 80357.5428571428) + end + + @test isapprox(energy3[], 0.0) + end + end +end diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 91ec13ef13..8b17226d88 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -41,9 +41,118 @@ @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end + + @trixi_testset "structure/oscillating_beam_2d.jl with EnergyCalculatorCallback" begin + # Load variables from the example + trixi_include(@__MODULE__, + joinpath(examples_dir(), "structure", "oscillating_beam_2d.jl"), + ode=nothing, sol=nothing) + + # We simply clamp all particles, move them up against gravity, and verify that + # the energy calculated is just the potential energy difference. + movement_function(x, t) = x + SVector(0.0, t) + is_moving(t) = true + prescribed_motion = PrescribedMotion(movement_function, is_moving) + + structure_system = TotalLagrangianSPHSystem(structure, smoothing_kernel, + smoothing_length, + material.E, material.nu, + n_clamped_particles=nparticles(structure), + acceleration=(0.0, -gravity), + clamped_particles_motion=prescribed_motion) + + semi = Semidiscretization(structure_system) + ode = semidiscretize(semi, (0.0, 1.0)) + + energy_calculator = EnergyCalculatorCallback{Float64}(structure_system, semi; + interval=1) + + sol = @trixi_test_nowarn solve(ode, RDPK3SpFSAL49(), save_everystep=false, + callback=energy_calculator) + + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + + # Potential energy difference should be m * g * h + @test isapprox(calculated_energy(energy_calculator), + sum(structure_system.mass) * gravity * 1) + end end @testset verbose=true "FSI" begin + @trixi_testset "fluid/hydrostatic_water_column_2d.jl with moving TLSPH walls" begin + # In this test, we move a water-filled tank up against gravity by 1 unit + # and verify that the energy calculated by the `EnergyCalculatorCallback` + # matches the expected potential energy. + + # Load variables from the example + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "hydrostatic_water_column_2d.jl"), + # Use tank without airspace + tank_size=(1.0, 0.9), + # Higher speed of sound for smaller error + # due to compressibility. + sound_speed=50.0, + ode=nothing, sol=nothing) + + # Move the tank up against gravity by 1 unit over 1 second + function movement_function(x, t) + return x + SVector(0.0, 0.5 * sin(pi * (t - 0.5)) + 0.5) + end + is_moving(t) = true + prescribed_motion = PrescribedMotion(movement_function, is_moving) + + # Create TLSPH system for the tank walls and clamp all particles. + # This is identical to a `WallBoundarySystem`, but now we can + # use the `EnergyCalculatorCallback` to compute the energy. + boundary_spacing = tank.boundary.particle_spacing + tlsph_kernel = WendlandC2Kernel{2}() + tlsph_smoothing_length = sqrt(2) * boundary_spacing + + tlsph_system = TotalLagrangianSPHSystem(tank.boundary, tlsph_kernel, + tlsph_smoothing_length, + 1.0e6, 0.3, + n_clamped_particles=nparticles(tank.boundary), + acceleration=(0.0, -gravity), + clamped_particles_motion=prescribed_motion, + boundary_model=boundary_model) + + semi = Semidiscretization(fluid_system, tlsph_system, + parallelization_backend=PolyesterBackend()) + ode = semidiscretize(semi, (0.0, 1.0)) + + # Energy calculators for fluid + tank and fluid only + energy_calculator1 = EnergyCalculatorCallback{eltype(tlsph_system)}(tlsph_system, + semi; + interval=1) + energy_calculator2 = EnergyCalculatorCallback{eltype(tlsph_system)}(tlsph_system, + semi; + interval=1, + only_compute_force_on_fluid=true) + + sol = @trixi_test_nowarn solve(ode, RDPK3SpFSAL35(), save_everystep=false, + callback=CallbackSet(info_callback, + energy_calculator1, + energy_calculator2)) + + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + + # Potential energy difference should be m * g * h and h = 1 + expected_energy_fluid = sum(fluid_system.mass) * gravity * 1 + expected_energy_tank = sum(tlsph_system.mass) * gravity * 1 + + # We don't expect very accurate results here because the fluid is weakly + # compressible and is deformed during the simulation. + # A slower prescribed motion (e.g., over 2 seconds instead of 1) or a higher + # speed of sound in the fluid would improve accuracy (and increase runtime). + @test isapprox(calculated_energy(energy_calculator1), + expected_energy_fluid + expected_energy_tank, rtol=5e-4) + @test isapprox(calculated_energy(energy_calculator2), expected_energy_fluid, + rtol=5e-4) + end + @trixi_testset "fsi/falling_water_column_2d.jl" begin @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "fsi",