Add energy calculator system

examples/fluid/oscillating_drop_2d.jl

@@ -7,6 +7,7 @@
 
 using TrixiParticles
 using OrdinaryDiffEq
+using LinearAlgebra: dot
 
 # ==========================================================================================
 # ==== Resolution
@@ -26,8 +27,8 @@ period = 4.567375
 tspan = (0.0, 1period)
 
 fluid_density = 1000.0
-sound_speed = 10.0
-state_equation = StateEquationCole(; sound_speed, exponent=7,
+sound_speed = 15 * sigma * radius
+state_equation = StateEquationCole(; sound_speed, exponent=1,
                                    reference_density=fluid_density)
 
 # Equation A.19 in the paper rearranged.
@@ -56,13 +57,81 @@ fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator,
                                            density_diffusion=density_diffusion,
                                            source_terms=source_terms)
 
+energy_system = TrixiParticles.EnergyCalculatorSystem{2}()
+
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system)
+semi = Semidiscretization(fluid_system, energy_system)
 ode = semidiscretize(semi, tspan)
 
+function kinetic_energy(v, u, t, system)
+    return nothing
+end
+
+function kinetic_energy(v, u, t, system::WeaklyCompressibleSPHSystem)
+    return sum(TrixiParticles.eachparticle(system)) do particle
+        vel = TrixiParticles.current_velocity(v, system, particle)
+        return system.mass[particle] * dot(vel, vel) / 2
+    end
+end
+
+function potential_energy(v, u, t, system)
+    return nothing
+end
+
+function potential_energy(v, u, t, system::WeaklyCompressibleSPHSystem)
+    return sum(TrixiParticles.eachparticle(system)) do particle
+        coords = TrixiParticles.current_coords(u, system, particle)
+        return system.mass[particle] * 0.5 * omega^2 * dot(coords, coords)
+    end
+end
+
+function mechanical_energy_normalized(v, u, t, system)
+    return nothing
+end
+
+function mechanical_energy_normalized(v, u, t, system::WeaklyCompressibleSPHSystem)
+    return (mechanical_energy(v, u, t, system) - 898.547) / 898.547
+end
+
+function mechanical_energy(v, u, t, system)
+    return nothing
+end
+
+function mechanical_energy(v, u, t, system::WeaklyCompressibleSPHSystem)
+    return kinetic_energy(v, u, t, system) + potential_energy(v, u, t, system)
+end
+
+function internal_energy(v, u, t, system)
+    return nothing
+end
+
+function internal_energy(v, u, t, system::WeaklyCompressibleSPHSystem)
+    return energy_system.current_energy[1]
+end
+
+function total_energy(v, u, t, system)
+    return nothing
+end
+
+function total_energy(v, u, t, system::WeaklyCompressibleSPHSystem)
+    return internal_energy(v, u, t, system) + mechanical_energy(v, u, t, system)
+end
+
+function total_energy_normalized(v, u, t, system)
+    return nothing
+end
+
+function total_energy_normalized(v, u, t, system::WeaklyCompressibleSPHSystem)
+    return (total_energy(v, u, t, system) - 898.547) / 898.547
+end
+
 info_callback = InfoCallback(interval=50)
-saving_callback = SolutionSavingCallback(dt=0.04, prefix="")
+saving_callback = SolutionSavingCallback(dt=0.1, prefix="";
+                                         kinetic_energy, potential_energy,
+                                         mechanical_energy, mechanical_energy_normalized,
+                                         internal_energy, total_energy,
+                                         total_energy_normalized)
 
 callbacks = CallbackSet(info_callback, saving_callback)
 

src/TrixiParticles.jl

@@ -34,7 +34,7 @@ using TrixiBase: trixi_include, @trixi_timeit, timer, timeit_debug_enabled,
                                 ParallelUpdate, SemiParallelUpdate, SerialUpdate
 using PointNeighbors: PointNeighbors, foreach_point_neighbor, copy_neighborhood_search,
                       @threaded
-using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, paraview_collection, vtk_save
+using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, paraview_collection, vtk_save, VTKFieldData
 
 # `util.jl` depends on the `GPUSystem` type defined in `system.jl`
 include("general/system.jl")

src/general/energy_calculator_system.jl

@@ -0,0 +1,93 @@
+struct EnergyCalculatorSystem{NDIMS, ELTYPE <: Real} <: System{NDIMS}
+    initial_energy::ELTYPE
+    current_energy::Vector{ELTYPE}
+
+    function EnergyCalculatorSystem{NDIMS}(; initial_energy=0.0) where {NDIMS}
+        new{NDIMS, typeof(initial_energy)}(initial_energy, [initial_energy])
+    end
+end
+
+timer_name(::EnergyCalculatorSystem) = "dummy"
+
+@inline Base.eltype(::EnergyCalculatorSystem{NDIMS, ELTYPE}) where {NDIMS, ELTYPE} = ELTYPE
+
+@inline nparticles(::EnergyCalculatorSystem) = 1
+
+@inline v_nvariables(::EnergyCalculatorSystem) = 1
+
+@inline initial_coordinates(::EnergyCalculatorSystem) = nothing
+
+function write_u0!(u0, ::EnergyCalculatorSystem)
+    # `u` can't be of size 0, or we have to write an extra `write2vtk` function
+    u0 .= 0
+
+    return u0
+end
+
+function write_v0!(v0, system::EnergyCalculatorSystem)
+    v0[1] = system.initial_energy
+
+    return v0
+end
+
+function nhs_coords(system, neighbor::EnergyCalculatorSystem, u)
+    return u
+end
+
+function update_quantities!(system::EnergyCalculatorSystem, v, u, v_ode, u_ode, semi, t)
+    system.current_energy[1] = v[1]
+
+    return system
+end
+
+@inline function interact!(dv_ode, v_ode, u_ode, system::EnergyCalculatorSystem,
+                           neighbor_system, semi; timer_str="")
+    @trixi_timeit timer() "calculate energy" begin
+        dv = wrap_v(dv_ode, system, semi)
+        dv_neighbor = wrap_v(dv_ode, neighbor_system, semi)
+        v_neighbor = wrap_v(v_ode, neighbor_system, semi)
+
+        for particle in eachparticle(neighbor_system)
+            d_velocity = current_velocity(dv_neighbor, neighbor_system, particle)
+            velocity = current_velocity(v_neighbor, neighbor_system, particle)
+
+            dv[1] -= dot(hydrodynamic_mass(neighbor_system, particle) * d_velocity,
+                         velocity)
+        end
+    end
+
+    return dv_ode
+end
+
+@inline function interact!(dv_ode, v_ode, u_ode, system::EnergyCalculatorSystem,
+                           neighbor_system::EnergyCalculatorSystem, semi; timer_str="")
+    return dv_ode
+end
+
+@inline function interact!(dv_ode, v_ode, u_ode, system,
+                           neighbor_system::EnergyCalculatorSystem, semi; timer_str="")
+    return dv_ode
+end
+
+vtkname(system::EnergyCalculatorSystem) = "energy"
+
+function write2vtk!(vtk, v, u, t, system::EnergyCalculatorSystem;
+                    write_meta_data=true)
+    vtk["energy", VTKFieldData()] = v[1]
+    vtk["initial_energy"] = system.initial_energy
+
+    return vtk
+end
+
+function Base.show(io::IO, system::EnergyCalculatorSystem)
+    print(io, "EnergyCalculatorSystem")
+end
+
+function Base.show(io::IO, ::MIME"text/plain", system::EnergyCalculatorSystem)
+    if get(io, :compact, false)
+        show(io, system)
+    else
+        summary_header(io, "EnergyCalculatorSystem")
+        summary_footer(io)
+    end
+end

src/general/general.jl

@@ -9,3 +9,4 @@ include("buffer.jl")
 include("interpolation.jl")
 include("custom_quantities.jl")
 include("neighborhood_search.jl")
+include("energy_calculator_system.jl")

src/general/semidiscretization.jl

@@ -189,6 +189,22 @@ end
     return compact_support(smoothing_kernel, smoothing_length)
 end
 
+@inline function compact_support(system::EnergyCalculatorSystem, neighbor)
+    # This NHS is never used
+    return 0.0
+end
+
+@inline function compact_support(system, neighbor::EnergyCalculatorSystem)
+    # This NHS is never used
+    return 0.0
+end
+
+@inline function compact_support(system::EnergyCalculatorSystem,
+                                 neighbor::EnergyCalculatorSystem)
+    # This NHS is never used
+    return 0.0
+end
+
 @inline function get_neighborhood_search(system, semi)
     (; neighborhood_searches) = semi
 
@@ -440,6 +456,7 @@ end
 end
 
 @inline add_velocity!(du, v, particle, system::BoundarySPHSystem) = du
+@inline add_velocity!(du, v, particle, system::EnergyCalculatorSystem) = du
 
 function kick!(dv_ode, v_ode, u_ode, semi, t)
     @trixi_timeit timer() "kick!" begin
@@ -817,6 +834,24 @@ function update!(neighborhood_search, system::GPUSystem, x, y; points_moving=(tr
                            parallelization_backend=KernelAbstractions.get_backend(system))
 end
 
+function nhs_coords(system::EnergyCalculatorSystem,
+                    neighbor, u)
+    # Don't update
+    return nothing
+end
+
+function nhs_coords(system,
+                    neighbor::EnergyCalculatorSystem, u)
+    # Don't update
+    return nothing
+end
+
+function nhs_coords(system::EnergyCalculatorSystem,
+                    neighbor::EnergyCalculatorSystem, u)
+    # Don't update
+    return nothing
+end
+
 function check_configuration(systems)
     foreach_system(systems) do system
         check_configuration(system, systems)