From 25b0d3e7f220d4858d13ff37ce19010db26ffa87 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 8 Nov 2024 00:03:02 +0100 Subject: [PATCH 1/6] Avoid bounds checking where it is safe to do so --- src/TrixiParticles.jl | 1 + src/general/density_calculators.jl | 6 +- src/general/system.jl | 24 ++++--- src/schemes/fluid/fluid.jl | 2 +- src/schemes/fluid/viscosity.jl | 36 ++++++---- .../fluid/weakly_compressible_sph/rhs.jl | 70 +++++++++++-------- .../fluid/weakly_compressible_sph/system.jl | 3 +- 7 files changed, 82 insertions(+), 60 deletions(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index bcdce98ac8..f5afcc5759 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -3,6 +3,7 @@ module TrixiParticles using Reexport: @reexport using Adapt: Adapt +using Base: @propagate_inbounds using CSV: CSV using Dates using DataFrames: DataFrame diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl index 74a5f616de..2a3d45ea09 100644 --- a/src/general/density_calculators.jl +++ b/src/general/density_calculators.jl @@ -23,15 +23,15 @@ difference of the coordinates, ``v_{ab} = v_a - v_b`` of the velocities of parti """ struct ContinuityDensity end -@inline function particle_density(v, system, particle) +@propagate_inbounds function particle_density(v, system, particle) particle_density(v, system.density_calculator, system, particle) end -@inline function particle_density(v, ::SummationDensity, system, particle) +@propagate_inbounds function particle_density(v, ::SummationDensity, system, particle) return system.cache.density[particle] end -@inline function particle_density(v, ::ContinuityDensity, system, particle) +@propagate_inbounds function particle_density(v, ::ContinuityDensity, system, particle) return v[end, particle] end diff --git a/src/general/system.jl b/src/general/system.jl index 674555aea4..9258672277 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -54,28 +54,31 @@ initialize!(system, neighborhood_search) = system @inline active_particles(system, ::Nothing) = eachparticle(system) # This should not be dispatched by system type. We always expect to get a column of `A`. -@inline function extract_svector(A, system, i) +@propagate_inbounds function extract_svector(A, system, i) extract_svector(A, Val(ndims(system)), i) end # Return the `i`-th column of the array `A` as an `SVector`. @inline function extract_svector(A, ::Val{NDIMS}, i) where {NDIMS} - return SVector(ntuple(@inline(dim->A[dim, i]), NDIMS)) + # Explicit bounds check, which can be removed by calling this function with `@inbounds` + @boundscheck checkbounds(A, NDIMS, i) + + # Assume inbounds access now + return SVector(ntuple(@inline(dim -> @inbounds A[dim, i]), NDIMS)) end # Return `A[:, :, i]` as an `SMatrix`. @inline function extract_smatrix(A, system, particle) # Extract the matrix elements for this particle as a tuple to pass to SMatrix - return SMatrix{ndims(system), ndims(system)}( - # Convert linear index to Cartesian index - ntuple(@inline(i->A[mod(i - 1, ndims(system)) + 1, - div(i - 1, ndims(system)) + 1, - particle]), - Val(ndims(system)^2))) + return SMatrix{ndims(system), + ndims(system)}(ntuple(@inline(i->A[mod(i - 1, ndims(system)) + 1, + div(i - 1, ndims(system)) + 1, + particle]), + Val(ndims(system)^2))) end # Specifically get the current coordinates of a particle for all system types. -@inline function current_coords(u, system, particle) +@propagate_inbounds function current_coords(u, system, particle) return extract_svector(current_coordinates(u, system), system, particle) end @@ -91,7 +94,8 @@ end # This can be dispatched by system type. @inline initial_coordinates(system) = system.initial_condition.coordinates -@inline current_velocity(v, system, particle) = extract_svector(v, system, particle) +@propagate_inbounds current_velocity(v, system, particle) = extract_svector(v, system, + particle) @inline function current_acceleration(system, particle) # TODO: Return `dv` of solid particles diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index e0f07796cb..fd0a9c3a56 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -13,7 +13,7 @@ function create_cache_density(ic, ::ContinuityDensity) return (;) end -@inline hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle] +@propagate_inbounds hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle] function write_u0!(u0, system::FluidSystem) (; initial_condition) = system diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 09bde77f9f..95c9ff88e8 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -1,9 +1,9 @@ # Unpack the neighboring systems viscosity to dispatch on the viscosity type -@inline function dv_viscosity(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) +@propagate_inbounds function dv_viscosity(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) viscosity = viscosity_model(particle_system, neighbor_system) return dv_viscosity(viscosity, particle_system, neighbor_system, @@ -12,10 +12,10 @@ sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) end -@inline function dv_viscosity(viscosity, particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) +@propagate_inbounds function dv_viscosity(viscosity, particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) return viscosity(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, @@ -105,12 +105,16 @@ function kinematic_viscosity(system, viscosity::ViscosityMorris) return viscosity.nu end -@inline function (viscosity::Union{ArtificialViscosityMonaghan, - ViscosityMorris})(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, - distance, sound_speed, m_a, m_b, - rho_a, rho_b, grad_kernel) +@propagate_inbounds function (viscosity::Union{ArtificialViscosityMonaghan, + ViscosityMorris})(particle_system, + neighbor_system, + v_particle_system, + v_neighbor_system, + particle, neighbor, + pos_diff, distance, + sound_speed, + m_a, m_b, rho_a, rho_b, + grad_kernel) (; smoothing_length) = particle_system rho_mean = 0.5 * (rho_a + rho_b) @@ -250,4 +254,6 @@ function kinematic_viscosity(system, viscosity::ViscosityAdami) return viscosity.nu end -@inline viscous_velocity(v, system, particle) = current_velocity(v, system, particle) +@propagate_inbounds function viscous_velocity(v, system, particle) + return current_velocity(v, system, particle) +end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 9e40eacfa1..c74b259510 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -23,19 +23,23 @@ function interact!(dv, v_particle_system, u_particle_system, foreach_point_neighbor(particle_system, neighbor_system, system_coords, neighbor_system_coords, neighborhood_search) do particle, neighbor, pos_diff, distance - rho_a = particle_density(v_particle_system, particle_system, particle) - rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + # `foreach_point_neighbor` makes sure that `particle` and `neighbor` are + # in bounds of the respective system. For performance reasons, we use `@inbounds` + # in this hot loop to avoid bounds checking when extracting particle quantities. + rho_a = @inbounds particle_density(v_particle_system, particle_system, particle) + rho_b = @inbounds particle_density(v_neighbor_system, neighbor_system, neighbor) rho_mean = 0.5 * (rho_a + rho_b) - # Determine correction values - viscosity_correction, pressure_correction, surface_tension_correction = free_surface_correction(correction, - particle_system, - rho_mean) + # Determine correction factors. + # This can be ignored, as these are all 1 when no correction is used. + (viscosity_correction, pressure_correction, + surface_tension_correction) = free_surface_correction(correction, particle_system, + rho_mean) grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) - m_a = hydrodynamic_mass(particle_system, particle) - m_b = hydrodynamic_mass(neighbor_system, neighbor) + m_a = @inbounds hydrodynamic_mass(particle_system, particle) + m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) # The following call is equivalent to # `p_a = particle_pressure(v_particle_system, particle_system, particle)` @@ -43,20 +47,23 @@ function interact!(dv, v_particle_system, u_particle_system, # Only when the neighbor system is a `BoundarySPHSystem` or a `TotalLagrangianSPHSystem` # with the boundary model `PressureMirroring`, this will return `p_b = p_a`, which is # the pressure of the fluid particle. - p_a, p_b = particle_neighbor_pressure(v_particle_system, v_neighbor_system, - particle_system, neighbor_system, - particle, neighbor) + p_a, p_b = @inbounds particle_neighbor_pressure(v_particle_system, + v_neighbor_system, + particle_system, neighbor_system, + particle, neighbor) dv_pressure = pressure_correction * pressure_acceleration(particle_system, neighbor_system, neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, distance, grad_kernel, correction) + # Propagate `@inbounds` to the viscosity function, which accesses particle data dv_viscosity_ = viscosity_correction * - dv_viscosity(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + @inbounds dv_viscosity(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, + grad_kernel) dv_surface_tension = surface_tension_correction * surface_tension_force(surface_tension_a, surface_tension_b, @@ -66,17 +73,19 @@ function interact!(dv, v_particle_system, u_particle_system, dv_adhesion = adhesion_force(surface_tension, particle_system, neighbor_system, particle, neighbor, pos_diff, distance) - @inbounds for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_surface_tension[i] + - dv_adhesion[i] + for i in 1:ndims(particle_system) + @inbounds dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + + dv_surface_tension[i] + dv_adhesion[i] # Debug example # debug_array[i, particle] += dv_pressure[i] end # TODO If variable smoothing_length is used, this should use the neighbor smoothing length - continuity_equation!(dv, density_calculator, v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, m_b, rho_a, rho_b, - particle_system, neighbor_system, grad_kernel) + # Propagate `@inbounds` to the continuity equation, which accesses particle data + @inbounds continuity_equation!(dv, density_calculator, v_particle_system, + v_neighbor_system, particle, neighbor, + pos_diff, distance, m_b, rho_a, rho_b, + particle_system, neighbor_system, grad_kernel) end # Debug example # periodic_box = neighborhood_search.periodic_box @@ -98,12 +107,12 @@ end end # This formulation was chosen to be consistent with the used pressure_acceleration formulations. -@inline function continuity_equation!(dv, density_calculator::ContinuityDensity, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, - particle_system::WeaklyCompressibleSPHSystem, - neighbor_system, grad_kernel) +@propagate_inbounds function continuity_equation!(dv, density_calculator::ContinuityDensity, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + m_b, rho_a, rho_b, + particle_system::WeaklyCompressibleSPHSystem, + neighbor_system, grad_kernel) (; density_diffusion) = particle_system vdiff = current_velocity(v_particle_system, particle_system, particle) - @@ -121,9 +130,10 @@ end end end -@inline function particle_neighbor_pressure(v_particle_system, v_neighbor_system, - particle_system, neighbor_system, - particle, neighbor) +@propagate_inbounds function particle_neighbor_pressure(v_particle_system, + v_neighbor_system, + particle_system, neighbor_system, + particle, neighbor) p_a = particle_pressure(v_particle_system, particle_system, particle) p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index eb6c06c849..6d66ba06a2 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -211,7 +211,8 @@ end return ndims(system) + 1 end -@inline function particle_pressure(v, system::WeaklyCompressibleSPHSystem, particle) +@propagate_inbounds function particle_pressure(v, system::WeaklyCompressibleSPHSystem, + particle) return system.pressure[particle] end From f31837709666b7e1275aa7ccff98a3cdecbc0d20 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 12 Nov 2024 17:49:35 +0100 Subject: [PATCH 2/6] Avoid one more bounds check in density diffusion --- src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index c4ddc8b819..7e4576420a 100644 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -199,7 +199,7 @@ function update!(density_diffusion::DensityDiffusionAntuono, neighborhood_search return density_diffusion end -@inline function density_diffusion!(dv, density_diffusion::DensityDiffusion, +@propagate_inbounds function density_diffusion!(dv, density_diffusion::DensityDiffusion, v_particle_system, particle, neighbor, pos_diff, distance, m_b, rho_a, rho_b, particle_system::FluidSystem, grad_kernel) From 22329dbf7c710255575c7a652b1d402030d1ee39 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 12 Nov 2024 18:01:30 +0100 Subject: [PATCH 3/6] Reformat code --- .../fluid/weakly_compressible_sph/density_diffusion.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index 7e4576420a..5d712235f8 100644 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -200,9 +200,9 @@ function update!(density_diffusion::DensityDiffusionAntuono, neighborhood_search end @propagate_inbounds function density_diffusion!(dv, density_diffusion::DensityDiffusion, - v_particle_system, particle, neighbor, pos_diff, - distance, m_b, rho_a, rho_b, - particle_system::FluidSystem, grad_kernel) + v_particle_system, particle, neighbor, + pos_diff, distance, m_b, rho_a, rho_b, + particle_system::FluidSystem, grad_kernel) # Density diffusion terms are all zero for distance zero distance < sqrt(eps()) && return From 09ab7bac2f8e3c95fc37fdf6809503d5362383e0 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 12 Nov 2024 18:49:07 +0100 Subject: [PATCH 4/6] Use `extract_svector` from PointNeighbors.jl --- src/general/system.jl | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/src/general/system.jl b/src/general/system.jl index 9258672277..cddde02100 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -55,16 +55,7 @@ initialize!(system, neighborhood_search) = system # This should not be dispatched by system type. We always expect to get a column of `A`. @propagate_inbounds function extract_svector(A, system, i) - extract_svector(A, Val(ndims(system)), i) -end - -# Return the `i`-th column of the array `A` as an `SVector`. -@inline function extract_svector(A, ::Val{NDIMS}, i) where {NDIMS} - # Explicit bounds check, which can be removed by calling this function with `@inbounds` - @boundscheck checkbounds(A, NDIMS, i) - - # Assume inbounds access now - return SVector(ntuple(@inline(dim -> @inbounds A[dim, i]), NDIMS)) + PointNeighbors.extract_svector(A, Val(ndims(system)), i) end # Return `A[:, :, i]` as an `SMatrix`. From 1fbb07aaab9c8688912a484cf8c70bf741e4ea7e Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 12 Nov 2024 18:55:25 +0100 Subject: [PATCH 5/6] Revert 09ab7bac2f8e3c95fc37fdf6809503d5362383e0 --- src/general/system.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/general/system.jl b/src/general/system.jl index cddde02100..48d556eba9 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -55,7 +55,16 @@ initialize!(system, neighborhood_search) = system # This should not be dispatched by system type. We always expect to get a column of `A`. @propagate_inbounds function extract_svector(A, system, i) - PointNeighbors.extract_svector(A, Val(ndims(system)), i) + extract_svector(A, Val(ndims(system)), i) +end + +# Return the `i`-th column of the array `A` as an `SVector`. +@inline function extract_svector(A, ::Val{NDIMS}, i) where {NDIMS} + # Explicit bounds check, which can be removed by calling this function with `@inbounds` + @boundscheck checkbounds(A, NDIMS, i) + + # Assume inbounds access now + return SVector(ntuple(@inline(dim->@inbounds A[dim, i]), NDIMS)) end # Return `A[:, :, i]` as an `SMatrix`. From 3adb4ea5055e58ae0e1328feefbd9b6d58134a84 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 12 Nov 2024 18:55:59 +0100 Subject: [PATCH 6/6] Fix tests --- test/systems/solid_system.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/systems/solid_system.jl b/test/systems/solid_system.jl index e5f0c0e057..53cf860dd3 100644 --- a/test/systems/solid_system.jl +++ b/test/systems/solid_system.jl @@ -131,7 +131,7 @@ Base.ntuple(f, ::Symbol) = ntuple(f, 2) # Make `extract_svector` work function TrixiParticles.current_coords(system::Val{:mock_system_tensor}, particle) - return TrixiParticles.extract_svector(current_coordinates[i], system, + return TrixiParticles.extract_svector(current_coordinates[i], Val(2), particle) end