Add EnergyCalculatorCallback

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,

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

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

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

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

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

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

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]

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