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
 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,207 @@
+
+"""
+    EnergyCalculatorCallback{ELTYPE}(; interval=1)
+
+Callback to calculate the energy contribution from clamped particles in
+[`TotalLagrangianSPHSystem`](@ref)s.
+
+This callback computes the work done by clamped particles that are moving according to a
+[`PrescribedMotion`](@ref).
+The energy is calculated by integrating the instantaneous power, which is the dot product
+of the force exerted by the clamped particle and its prescribed velocity.
+The callback uses a simple explicit Euler time integration for the energy calculation.
+
+!!! info "TLSPH System Configuration"
+    The [`TotalLagrangianSPHSystem`](@ref) must be created with
+    `use_with_energy_calculator_callback=true` to prepare the system for energy
+    calculation by allocating the necessary cache for clamped particles.
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
+
+# Arguments
+- `ELTYPE`: The floating point type used for time and energy calculations.
+            This should match the floating point type used in the
+            [`TotalLagrangianSPHSystem`](@ref), which can be obtained via `eltype(system)`.
+
+# Keywords
+- `interval=1`: Interval (in number of time steps) at which to compute the energy.
+                It is recommended to set this either to `1` (every time step) or to a small
+                integer (e.g., `2` or `5`) to reduce time integration errors
+                in the energy calculation.
+
+# Examples
+```jldoctest; output = false, setup = :(structure = RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); smoothing_kernel = WendlandC2Kernel{2}(); smoothing_length = 1.0; young_modulus = 1e6; poisson_ratio = 0.3; n_clamped_particles = 0; clamped_particles_motion = nothing)
+# Create TLSPH system with `use_with_energy_calculator_callback=true`
+system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothing_length,
+                                  young_modulus, poisson_ratio,
+                                  n_clamped_particles=n_clamped_particles,
+                                  clamped_particles_motion=clamped_particles_motion,
+                                  use_with_energy_calculator_callback=true) # Important!
+
+# Create an energy calculator that runs every 2 time steps
+energy_cb = EnergyCalculatorCallback{eltype(system)}(; interval=2)
+
+# output
+┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+│ EnergyCalculatorCallback{Float64}                                                                │
+│ ═════════════════════════════════                                                                │
+│ interval: ……………………………………………………… 2                                                                │
+└──────────────────────────────────────────────────────────────────────────────────────────────────┘
+```
+"""
+struct EnergyCalculatorCallback{T}
+    interval :: Int
+    t        :: T
+    energy   :: T
+end
+
+function EnergyCalculatorCallback{ELTYPE}(; interval=1) where {ELTYPE <: Real}
+    cb = EnergyCalculatorCallback(interval, Ref(zero(ELTYPE)), Ref(zero(ELTYPE)))
+
+    # The first one is the `condition`, the second the `affect!`
+    return DiscreteCallback(cb, cb, save_positions=(false, false))
+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
+
+    update_energy_calculator!(energy, v_ode, u_ode, semi.systems[1], semi, t, dt)
+
+    # Tell OrdinaryDiffEq that `u` has not been modified
+    u_modified!(integrator, false)
+
+    return integrator
+end
+
+function update_energy_calculator!(energy, v_ode, u_ode,
+                                   system::AbstractStructureSystem, 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
+
+        dv_clamped = system.cache.dv_clamped
+        set_zero!(dv_clamped)
+
+        dv = zero(wrap_v(v_ode, system, semi))
+        dv_combined = CombinedMatrix(dv, dv_clamped)
+        v = wrap_v(v_ode, system, semi)
+        u = wrap_u(u_ode, system, semi)
+        eachparticle = (n_integrated_particles(system) + 1):nparticles(system)
+
+        foreach_system(semi) do neighbor_system
+            v_neighbor = wrap_v(v_ode, neighbor_system, semi)
+            u_neighbor = wrap_u(u_ode, neighbor_system, semi)
+
+            interact!(dv_combined, v, u, v_neighbor, u_neighbor,
+                      system, neighbor_system, semi,
+                      eachparticle=eachparticle)
+        end
+
+        @threaded semi for particle in eachparticle
+            add_acceleration!(dv_combined, particle, system)
+            add_source_terms_inner!(dv_combined, v, u, particle, system,
+                                    source_terms(system), t)
+        end
+
+        n_clamped_particles = nparticles(system) - n_integrated_particles(system)
+        for clamped_particle in 1:n_clamped_particles
+            particle = clamped_particle + n_integrated_particles(system)
+            velocity = current_velocity(nothing, system, particle)
+            dv_particle = extract_svector(dv_clamped, system, clamped_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 done BY the particle times prescribed velocity.
+            # The work done BY the particle is the negative of the work done ON it.
+            energy[] -= dot(F_particle, velocity) * dt
+        end
+    end
+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{T}}) where {T}
+    @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
+
+# Data type that combines two matrices to behave like a single matrix.
+# This is used above to combine the `dv` for integrated particles and the `dv_clamped`
+# for clamped particles into a single matrix that can be passed to `interact!`.
+struct CombinedMatrix{T, M1, M2} <: AbstractMatrix{T}
+    matrix1::M1
+    matrix2::M2
+
+    function CombinedMatrix(matrix1, matrix2)
+        @assert size(matrix1, 1) == size(matrix2, 1)
+        @assert eltype(matrix1) == eltype(matrix2)
+
+        new{eltype(matrix1), typeof(matrix1), typeof(matrix2)}(matrix1, matrix2)
+    end
+end
+
+@inline function Base.size(cm::CombinedMatrix)
+    return (size(cm.matrix1, 1), size(cm.matrix1, 2) + size(cm.matrix2, 2))
+end
+
+@inline function Base.getindex(cm::CombinedMatrix, i, j)
+    @boundscheck checkbounds(cm, i, j)
+
+    length1 = size(cm.matrix1, 2)
+    if j <= length1
+        return @inbounds cm.matrix1[i, j]
+    else
+        return @inbounds cm.matrix2[i, j - length1]
+    end
+end
+
+@inline function Base.setindex!(cm::CombinedMatrix, value, i, j)
+    @boundscheck checkbounds(cm, i, j)
+
+    length1 = size(cm.matrix1, 2)
+    if j <= length1
+        return @inbounds cm.matrix1[i, j] = value
+    else
+        return @inbounds cm.matrix2[i, j - length1] = value
+    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/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

@@ -91,6 +91,7 @@ end
 function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length,
                                   young_modulus, poisson_ratio;
                                   n_clamped_particles=0, clamped_particles_motion=nothing,
+                                  use_with_energy_calculator_callback=false,
                                   boundary_model=nothing,
                                   acceleration=ntuple(_ -> 0.0,
                                                       ndims(smoothing_kernel)),
@@ -128,7 +129,9 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing
     initialize_prescribed_motion!(clamped_particles_motion, initial_condition,
                                   n_clamped_particles)
 
-    cache = create_cache_tlsph(clamped_particles_motion, initial_condition)
+    cache = (; create_cache_tlsph(clamped_particles_motion, initial_condition)...,
+             create_cache_tlsph(Val(use_with_energy_calculator_callback),
+                                initial_condition, n_clamped_particles)...)
 
     return TotalLagrangianSPHSystem(initial_condition, initial_coordinates,
                                     current_coordinates, mass, correction_matrix,
@@ -149,6 +152,15 @@ function create_cache_tlsph(::PrescribedMotion, initial_condition)
     return (; velocity, acceleration)
 end
 
+create_cache_tlsph(::Val{false}, initial_condition, n_clamped_particles) = (;)
+
+function create_cache_tlsph(::Val{true}, initial_condition, n_clamped_particles)
+    dv_clamped = Array{eltype(initial_condition)}(undef, ndims(initial_condition),
+                                                  n_clamped_particles)
+
+    return (; dv_clamped)
+end
+
 @inline function Base.eltype(::TotalLagrangianSPHSystem{<:Any, <:Any, ELTYPE}) where {ELTYPE}
     return ELTYPE
 end
@@ -336,11 +348,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