diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl index 4e61f8c574..2353390575 100644 --- a/src/callbacks/split_integration.jl +++ b/src/callbacks/split_integration.jl @@ -203,9 +203,8 @@ function kick_split!(dv_ode_split, v_ode_split, u_ode_split, p, t) v_ode_split, u_ode_split, semi_split) end - @trixi_timeit timer() "source terms" add_source_terms!(dv_ode_split, v_ode_split, - u_ode_split, semi, t; - semi_wrap=semi_split) + add_source_terms!(dv_ode_split, v_ode_split, u_ode_split, semi, t; + semi_wrap=semi_split) end function drift_split!(du_ode, v_ode, u_ode, p, t) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index cd407858e4..c6c6568af1 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -413,62 +413,69 @@ end function drift!(du_ode, v_ode, u_ode, semi, t) @trixi_timeit timer() "drift!" begin - @trixi_timeit timer() "reset ∂u/∂t" set_zero!(du_ode) - - @trixi_timeit timer() "velocity" begin - # Set velocity and add acceleration for each system - foreach_system(semi) do system - du = wrap_u(du_ode, system, semi) - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - - integrate_tlsph = semi.integrate_tlsph[] - @threaded semi for particle in each_integrated_particle(system) - # This can be dispatched per system - add_velocity!(du, v, u, particle, system, integrate_tlsph, t) - end - end + foreach_system(semi) do system + du = wrap_u(du_ode, system, semi) + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + set_velocity!(du, v, u, system, semi, t) end end return du_ode end -@inline function add_velocity!(du, v, u, particle, system, integrate_tlsph, t) - add_velocity!(du, v, u, particle, system, t) -end - -@inline function add_velocity!(du, v, u, particle, system::TotalLagrangianSPHSystem, - integrate_tlsph, t) - # Only add velocity for TLSPH systems if they are integrated - if integrate_tlsph - add_velocity!(du, v, u, particle, system, t) - end +# Generic fallback for all systems that don't define this function +function set_velocity!(du, v, u, system, semi, t) + set_velocity_default!(du, v, u, system, semi, t) end -@inline function add_velocity!(du, v, u, particle, system, t) - # Generic fallback for all systems that don't define this function - for i in 1:ndims(system) - @inbounds du[i, particle] = v[i, particle] +# Only set velocity for TLSPH systems if they are integrated +function set_velocity!(du, v, u, system::TotalLagrangianSPHSystem, semi, t) + if semi.integrate_tlsph[] + set_velocity_default!(du, v, u, system, semi, t) + else + set_zero!(du) end return du end # Solid wall boundary system doesn't integrate the particle positions -@inline add_velocity!(du, v, u, particle, system::WallBoundarySystem, t) = du +function set_velocity!(du, v, u, system::WallBoundarySystem, semi, t) + # Note that `du` is of length zero, so we don't have to set it to zero + return du +end -@inline function add_velocity!(du, v, u, particle, system::AbstractFluidSystem, t) - # This is zero unless a shifting technique is used - delta_v_ = delta_v(system, particle) +# Fluid systems integrate the particle positions and can have a shifting velocity +function set_velocity!(du, v, u, system::AbstractFluidSystem, semi, t) + @threaded semi for particle in each_integrated_particle(system) + delta_v_ = @inbounds delta_v(system, particle) - for i in 1:ndims(system) - @inbounds du[i, particle] = v[i, particle] + delta_v_[i] + for i in 1:ndims(system) + @inbounds du[i, particle] = v[i, particle] + delta_v_[i] + end end return du end +function set_velocity_default!(du, v, u, system, semi, t) + @threaded semi for particle in each_integrated_particle(system) + for i in 1:ndims(system) + @inbounds du[i, particle] = v[i, particle] + end + end + + return du +end + +# This defaults to optimized GPU copy that is about 4x faster than the threaded version above +function set_velocity_default!(du::AbstractGPUArray, v, u, system, semi, t) + indices = CartesianIndices(du) + copyto!(du, indices, v, indices) +end + function kick!(dv_ode, v_ode, u_ode, semi, t) @trixi_timeit timer() "kick!" begin # Check that the `UpdateCallback` is used if required @@ -482,8 +489,7 @@ function kick!(dv_ode, v_ode, u_ode, semi, t) @trixi_timeit timer() "system interaction" system_interaction!(dv_ode, v_ode, u_ode, semi) - @trixi_timeit timer() "source terms" add_source_terms!(dv_ode, v_ode, u_ode, - semi, t) + add_source_terms!(dv_ode, v_ode, u_ode, semi, t) end return dv_ode @@ -545,6 +551,7 @@ end # The `SplitIntegrationCallback` overwrites `semi_wrap` to use a different # semidiscretization for wrapping arrays. +# `semi_wrap` is the small semidiscretization, `semi` is the large semidiscretization. # TODO `semi` is not used yet, but will be used when the source terms API is modified # to match the custom quantities API. function add_source_terms!(dv_ode, v_ode, u_ode, semi, t; semi_wrap=semi) @@ -555,54 +562,63 @@ function add_source_terms!(dv_ode, v_ode, u_ode, semi, t; semi_wrap=semi) # `integrate_tlsph` is extracted from the `semi_wrap`, so that this function # can be used in the `SplitIntegrationCallback` as well. - integrate_tlsph = semi_wrap.integrate_tlsph[] - - @threaded semi for particle in each_integrated_particle(system) - # Dispatch by system type to exclude boundary systems - add_acceleration!(dv, particle, system, integrate_tlsph) - add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t, - integrate_tlsph) - end + # In this case, `semi_wrap` will be the small sub-integration semidiscretization. + add_source_terms!(dv, v, u, system, semi, t, semi_wrap.integrate_tlsph[]) end return dv_ode end -@inline source_terms(system) = nothing -@inline source_terms(system::Union{AbstractFluidSystem, AbstractStructureSystem}) = system.source_terms - -@inline function add_acceleration!(dv, particle, system, integrate_tlsph) - add_acceleration!(dv, particle, system) +# This is a no-op by default but can be dispatched by system type +function add_source_terms!(dv, v, u, system, semi, t, integrate_tlsph) + return dv end -@inline function add_acceleration!(dv, particle, system::TotalLagrangianSPHSystem, - integrate_tlsph) - integrate_tlsph && add_acceleration!(dv, particle, system) +function add_source_terms!(dv, v, u, + system::Union{AbstractFluidSystem, AbstractStructureSystem}, + semi, t, integrate_tlsph) + add_source_terms_inner!(dv, v, u, system, semi, t) end -@inline add_acceleration!(dv, particle, system) = dv +function add_source_terms!(dv, v, u, system::TotalLagrangianSPHSystem, + semi, t, integrate_tlsph) + if integrate_tlsph + add_source_terms_inner!(dv, v, u, system, semi, t) + end -@propagate_inbounds function add_acceleration!(dv, particle, - system::Union{AbstractFluidSystem, - AbstractStructureSystem}) - (; acceleration) = system + return dv +end - for i in 1:ndims(system) - dv[i, particle] += acceleration[i] +function add_source_terms_inner!(dv, v, u, + system::Union{AbstractFluidSystem, + AbstractStructureSystem}, + semi, t) + if iszero(system.acceleration) && isnothing(source_terms(system)) + # Nothing to do + return dv + end + + @trixi_timeit timer() "source terms" begin + @threaded semi for particle in each_integrated_particle(system) + add_acceleration!(dv, system, particle) + add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t) + end end return dv end -@inline function add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t, - integrate_tlsph) - add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t) -end +@inline source_terms(system) = nothing +@inline source_terms(system::Union{AbstractFluidSystem, AbstractStructureSystem}) = system.source_terms -@inline function add_source_terms_inner!(dv, v, u, particle, - system::TotalLagrangianSPHSystem, - source_terms_, t, integrate_tlsph) - integrate_tlsph && add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t) +@inline function add_acceleration!(dv, system, particle) + (; acceleration) = system + + for i in 1:ndims(system) + @inbounds dv[i, particle] += acceleration[i] + end + + return dv end @propagate_inbounds function add_source_terms_inner!(dv, v, u, particle, diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index d6439741ac..39ae4e9635 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -259,7 +259,13 @@ function kinetic_energy(system::ParticlePackingSystem, v_ode, u_ode, semi, t) end @inline source_terms(system::ParticlePackingSystem) = nothing -@inline add_acceleration!(dv, particle, system::ParticlePackingSystem) = dv + +# Special case for `ParticlePackingSystem`, which is an `AbstractFluidSystem` +# but doesn't have source terms or gravity/acceleration. +function add_source_terms!(dv, v, u, system::ParticlePackingSystem, + semi, t, integrate_tlsph) + return dv +end # Update from `UpdateCallback` (between time steps) update_particle_packing(system, v_ode, u_ode, semi, integrator) = system @@ -390,15 +396,19 @@ end end # Skip for fixed systems -@inline function add_velocity!(du, v, u, particle, - system::ParticlePackingSystem{<:Any, true}, t) +@inline function set_velocity!(du, v, u, + system::ParticlePackingSystem{<:Any, true}, + semi, t) + # Note that `du` is of size `(0, n_particles)`, so we don't have to set it to zero return du end # Add advection velocity -@inline function add_velocity!(du, v, u, particle, system::ParticlePackingSystem, t) - for i in 1:ndims(system) - du[i, particle] = system.advection_velocity[i, particle] +@inline function set_velocity!(du, v, u, system::ParticlePackingSystem, semi, t) + @threaded semi for particle in each_integrated_particle(system) + for i in 1:ndims(system) + @inbounds du[i, particle] = system.advection_velocity[i, particle] + end end return du diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 94aee1b048..427dd1632e 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -286,17 +286,20 @@ function calculate_dt(v_ode, u_ode, cfl_number, system::OpenBoundarySystem, semi return Inf end -@inline function add_velocity!(du, v, u, particle, system::OpenBoundarySystem, t) - boundary_zone = current_boundary_zone(system, particle) +@inline function set_velocity!(du, v, u, system::OpenBoundarySystem, semi, t) + @threaded semi for particle in each_integrated_particle(system) + boundary_zone = current_boundary_zone(system, particle) - pos = current_coords(u, system, particle) - v_particle = reference_velocity(boundary_zone, v, system, particle, pos, t) + pos = @inbounds current_coords(u, system, particle) + v_particle = @inbounds reference_velocity(boundary_zone, v, system, + particle, pos, t) - # This is zero unless a shifting technique is used - delta_v_ = delta_v(system, particle) + # This is zero unless a shifting technique is used + delta_v_ = @inbounds delta_v(system, particle) - for i in 1:ndims(system) - @inbounds du[i, particle] = v_particle[i] + delta_v_[i] + for i in 1:ndims(system) + @inbounds du[i, particle] = v_particle[i] + delta_v_[i] + end end return du diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl index ceaff12a21..a0760614e8 100644 --- a/src/schemes/fluid/shifting_techniques.jl +++ b/src/schemes/fluid/shifting_techniques.jl @@ -403,7 +403,9 @@ function update_shifting_from_callback!(system, v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) - update_shifting_inner!(system, shifting, v, u, v_ode, u_ode, semi) + @trixi_timeit timer() "update shifting" begin + update_shifting_inner!(system, shifting, v, u, v_ode, u_ode, semi) + end end # `ParticleShiftingTechnique{<:Any, <:Any, <:Any, true}` diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index cac4ed1f12..1e7b22cb92 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -419,6 +419,17 @@ function update_tlsph_positions!(system::TotalLagrangianSPHSystem, u, semi) return system end +# This defaults to optimized GPU copy that is about 4x faster than the threaded version above +function update_tlsph_positions!(system::TotalLagrangianSPHSystem, + u::AbstractGPUArray, semi) + (; current_coordinates) = system + + indices = CartesianIndices(u) + copyto!(current_coordinates, indices, u, indices) + + return system +end + function apply_prescribed_motion!(system::TotalLagrangianSPHSystem, prescribed_motion::PrescribedMotion, semi, t) (; clamped_particles_moving, current_coordinates, cache) = system diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index 3c2dcbeb58..050581687c 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -1,8 +1,8 @@ # Use `@trixi_testset` to isolate the mock functions in a separate namespace @trixi_testset "Semidiscretization" begin # Mock systems - struct System1 <: TrixiParticles.AbstractSystem{3} end - struct System2 <: TrixiParticles.AbstractSystem{3} end + struct System1 <: TrixiParticles.AbstractStructureSystem{3} end + struct System2 <: TrixiParticles.AbstractStructureSystem{3} end system1 = System1() system2 = System2() @@ -10,9 +10,9 @@ Base.eltype(::System1) = Float64 TrixiParticles.coordinates_eltype(::System1) = Float32 TrixiParticles.u_nvariables(::System1) = 3 - TrixiParticles.u_nvariables(::System2) = 4 + TrixiParticles.u_nvariables(::System2) = 3 TrixiParticles.v_nvariables(::System1) = 3 - TrixiParticles.v_nvariables(::System2) = 2 + TrixiParticles.v_nvariables(::System2) = 4 TrixiParticles.nparticles(::System1) = 2 TrixiParticles.nparticles(::System2) = 3 TrixiParticles.n_integrated_particles(::System1) = 2 @@ -25,8 +25,8 @@ semi = Semidiscretization(system1, system2, neighborhood_search=nothing) # Verification - @test semi.ranges_u == (1:6, 7:18) - @test semi.ranges_v == (1:6, 7:12) + @test semi.ranges_u == (1:6, 7:15) + @test semi.ranges_v == (1:6, 7:18) nhs = [TrixiParticles.TrivialNeighborhoodSearch{3}(search_radius=0.2, eachpoint=1:2) @@ -146,14 +146,27 @@ end @testset verbose=true "Source Terms" begin + function Base.getproperty(::System1, name::Symbol) + if name == :acceleration + return SVector(0.0, 0.0, 0.0) + end + error("property $name not defined for `System1`") + end + function Base.getproperty(::System2, name::Symbol) + if name == :acceleration + return SVector(0.0, 0.0, 1.0) + end + error("property $name not defined for `System2`") + end TrixiParticles.source_terms(::System1) = SourceTermDamping(damping_coefficient=0.1) + TrixiParticles.source_terms(::System2) = nothing TrixiParticles.current_density(v, system::System1, particle) = 0.0 TrixiParticles.current_pressure(v, system::System1, particle) = 0.0 semi = Semidiscretization(system1, system2, neighborhood_search=nothing) - dv_ode = zeros(3 * 2 + 2 * 3) - du_ode = zeros(3 * 2 + 4 * 3) + dv_ode = zeros(3 * 2 + 4 * 3) + du_ode = zeros(3 * 2 + 3 * 3) u_ode = zero(du_ode) v1 = [1.0 2.0 @@ -168,6 +181,31 @@ @test dv1 == -0.1 * v1 dv2 = TrixiParticles.wrap_v(dv_ode, system2, semi) - @test iszero(dv2) + @test dv2 == vcat(zeros(2, 3), ones(1, 3), zeros(1, 3)) + end + + @testset verbose=true "drift!" begin + semi = Semidiscretization(system1, system2, neighborhood_search=nothing) + + du_ode = fill(NaN, 3 * 2 + 3 * 3) + u_ode = zeros(3 * 2 + 3 * 3) + + v1 = [1.0 2.0 + 3.0 4.0 + 5.0 6.0] + v2 = [10.0 11.0 12.0 + 20.0 21.0 22.0 + 30.0 31.0 32.0 + -999.0 -999.0 -999.0] + v_ode = vcat(vec(v1), vec(v2)) + + returned = TrixiParticles.drift!(du_ode, v_ode, u_ode, semi, 0.0) + @test returned === du_ode + + du1 = TrixiParticles.wrap_u(du_ode, system1, semi) + @test du1 == v1 + + du2 = TrixiParticles.wrap_u(du_ode, system2, semi) + @test du2 == v2[1:3, :] end end diff --git a/test/schemes/structure/total_lagrangian_sph/rhs.jl b/test/schemes/structure/total_lagrangian_sph/rhs.jl index 69f4a920c4..6383f6977e 100644 --- a/test/schemes/structure/total_lagrangian_sph/rhs.jl +++ b/test/schemes/structure/total_lagrangian_sph/rhs.jl @@ -93,9 +93,6 @@ TrixiParticles.each_integrated_particle(::MockSystem) = each_integrated_particle TrixiParticles.smoothing_length(::MockSystem, _) = eps() - function TrixiParticles.add_acceleration!(_, _, ::MockSystem) - return nothing - end TrixiParticles.kernel_deriv(::Val{:mock_smoothing_kernel}, _, _) = kernel_deriv TrixiParticles.compact_support(::MockSystem, ::MockSystem) = 1000.0 function TrixiParticles.current_coords(system::MockSystem, particle)