From 05e4870988627946602d0aed31ec17dfe1b49ae8 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 21 May 2025 16:23:03 +0200 Subject: [PATCH 1/2] Make PST GPU-compatible --- src/callbacks/particle_shifting.jl | 16 ++++++++++------ .../fluid/weakly_compressible_sph/system.jl | 17 +++++++++++++++++ 2 files changed, 27 insertions(+), 6 deletions(-) diff --git a/src/callbacks/particle_shifting.jl b/src/callbacks/particle_shifting.jl index 145c9c1a1..c5d5a2156 100644 --- a/src/callbacks/particle_shifting.jl +++ b/src/callbacks/particle_shifting.jl @@ -32,7 +32,7 @@ function particle_shifting!(integrator) v_ode, u_ode = integrator.u.x dt = integrator.dt # Internal cache vector, which is safe to use as temporary array - u_cache = first(get_tmp_cache(integrator)) + vu_cache = first(get_tmp_cache(integrator)) # 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. @@ -41,7 +41,7 @@ function particle_shifting!(integrator) @trixi_timeit timer() "particle shifting" foreach_system(semi) do system u = wrap_u(u_ode, system, semi) v = wrap_v(v_ode, system, semi) - particle_shifting!(u, v, system, v_ode, u_ode, semi, u_cache, dt) + particle_shifting!(u, v, system, v_ode, u_ode, semi, vu_cache, dt) end # Tell OrdinaryDiffEq that `u` has been modified @@ -55,14 +55,18 @@ function particle_shifting!(u, v, system, v_ode, u_ode, semi, u_cache, dt) end function particle_shifting!(u, v, system::FluidSystem, v_ode, u_ode, semi, - u_cache, dt) + vu_cache, dt) # Wrap the cache vector to an NDIMS x NPARTICLES matrix. # We need this buffer because we cannot safely update `u` while iterating over it. + _, u_cache = vu_cache.x delta_r = wrap_u(u_cache, system, semi) set_zero!(delta_r) - v_max = maximum(particle -> norm(current_velocity(v, system, particle)), - eachparticle(system)) + # This has similar performance to `maximum(..., eachparticle(system))`, + # but is GPU-compatible. + v_max = maximum(x -> sqrt(dot(x, x)), + reinterpret(reshape, SVector{ndims(system), eltype(v)}, + current_velocity(v, system))) # TODO this needs to be adapted to multi-resolution. # Section 3.2 explains what else needs to be changed. @@ -88,7 +92,7 @@ function particle_shifting!(u, v, system::FluidSystem, v_ode, u_ode, semi, grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) # According to p. 29 below Eq. 9 - R = 0.2 + R = 2 // 10 n = 4 # Eq. 7 in Sun et al. (2017). diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 74306d844..aa8651dae 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -247,6 +247,23 @@ end system_correction(system::WeaklyCompressibleSPHSystem) = system.correction +@inline function current_velocity(v, system::WeaklyCompressibleSPHSystem) + return current_velocity(v, system.density_calculator, system) +end + +@inline function current_velocity(v, ::SummationDensity, + system::WeaklyCompressibleSPHSystem) + # When using `SummationDensity`, `v` contains only the velocity + return v +end + +@inline function current_velocity(v, ::ContinuityDensity, + system::WeaklyCompressibleSPHSystem) + # When using `ContinuityDensity`, the velocity is stored + # in the first `ndims(system)` rows of `v`. + return view(v, 1:ndims(system), :) +end + @inline function current_density(v, system::WeaklyCompressibleSPHSystem) return current_density(v, system.density_calculator, system) end From 689a6b951ac36dd2984a3b06a4a7ed3909cb3df8 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 21 May 2025 16:27:11 +0200 Subject: [PATCH 2/2] Reformat --- src/schemes/fluid/weakly_compressible_sph/system.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index aa8651dae..3cbd4b31a 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -252,13 +252,13 @@ system_correction(system::WeaklyCompressibleSPHSystem) = system.correction end @inline function current_velocity(v, ::SummationDensity, - system::WeaklyCompressibleSPHSystem) + system::WeaklyCompressibleSPHSystem) # When using `SummationDensity`, `v` contains only the velocity return v end @inline function current_velocity(v, ::ContinuityDensity, - system::WeaklyCompressibleSPHSystem) + system::WeaklyCompressibleSPHSystem) # When using `ContinuityDensity`, the velocity is stored # in the first `ndims(system)` rows of `v`. return view(v, 1:ndims(system), :)