Skip to content

Commit 995f182

Browse files
LasNikasLasNikasefaulhaber
authored
Make custom_quantities GPU compatible (#867)
* make it compatible * apply formatter * add tests * implement suggestions --------- Co-authored-by: LasNikas <niklas.nehe@web.de> Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
1 parent c9bb1ae commit 995f182

5 files changed

Lines changed: 169 additions & 9 deletions

File tree

‎src/callbacks/steady_state_reached.jl‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -26,17 +26,19 @@ end
2626

2727
function SteadyStateReachedCallback(; interval::Integer=0, dt=0.0,
2828
interval_size::Integer=10, abstol=1.0e-8, reltol=1.0e-6)
29+
ELTYPE = eltype(abstol)
2930
abstol, reltol = promote(abstol, reltol)
3031

3132
if dt > 0 && interval > 0
3233
throw(ArgumentError("setting both `interval` and `dt` is not supported"))
3334
end
3435

3536
if dt > 0
36-
interval = Float64(dt)
37+
interval = convert(ELTYPE, dt)
3738
end
3839

39-
steady_state_callback = SteadyStateReachedCallback(interval, abstol, reltol, [Inf64],
40+
steady_state_callback = SteadyStateReachedCallback(interval, abstol, reltol,
41+
[convert(ELTYPE, Inf)],
4042
interval_size)
4143

4244
if dt > 0

‎src/general/custom_quantities.jl‎

Lines changed: 11 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -6,22 +6,26 @@ Returns the total kinetic energy of all particles in a system.
66
function kinetic_energy(system, v_ode, u_ode, semi, t)
77
v = wrap_v(v_ode, system, semi)
88

9-
# If `each_moving_particle` is empty (no moving particles), return zero
10-
return sum(each_moving_particle(system), init=zero(eltype(system))) do particle
11-
velocity = current_velocity(v, system, particle)
12-
return 0.5 * system.mass[particle] * dot(velocity, velocity)
9+
# TODO: `current_velocity` should only contain active particles
10+
# (see https://github.com/trixi-framework/TrixiParticles.jl/issues/850)
11+
velocity = reinterpret(reshape, SVector{ndims(system), eltype(v)},
12+
view(current_velocity(v, system), :, active_particles(system)))
13+
mass = view(system.mass, active_particles(system))
14+
15+
return mapreduce(+, velocity, mass) do v_i, m_i
16+
return m_i * dot(v_i, v_i) / 2
1317
end
1418
end
1519

20+
kinetic_energy(system::BoundarySystem, v_ode, u_ode, semi, t) = zero(eltype(system))
21+
1622
"""
1723
total_mass
1824
1925
Returns the total mass of all particles in a system.
2026
"""
2127
function total_mass(system, v_ode, u_ode, semi, t)
22-
return sum(eachparticle(system)) do particle
23-
return system.mass[particle]
24-
end
28+
return sum(system.mass)
2529
end
2630

2731
function total_mass(system::BoundarySystem, v_ode, u_ode, semi, t)

‎test/examples/gpu.jl‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -354,6 +354,44 @@ end
354354
backend = TrixiParticles.KernelAbstractions.get_backend(sol.u[end].x[1])
355355
@test backend == Main.parallelization_backend
356356
end
357+
358+
@trixi_testset "fluid/pipe_flow_2d.jl - steady state reached (`dt`)" begin
359+
steady_state_reached = SteadyStateReachedCallback(; dt=0.002f0,
360+
interval_size=10,
361+
reltol=1.0f-3)
362+
363+
@trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__,
364+
joinpath(examples_dir(),
365+
"fluid",
366+
"pipe_flow_2d.jl"),
367+
extra_callback=steady_state_reached,
368+
tspan=(0.0f0, 1.5f0),
369+
parallelization_backend=Main.parallelization_backend,
370+
viscosity_boundary=nothing)
371+
372+
# Make sure that the simulation is terminated after a reasonable amount of time
373+
@test 0.1 < sol.t[end] < 1.0
374+
@test sol.retcode == ReturnCode.Terminated
375+
end
376+
377+
@trixi_testset "fluid/pipe_flow_2d.jl - steady state reached (`interval`)" begin
378+
steady_state_reached = SteadyStateReachedCallback(; interval=1,
379+
interval_size=10,
380+
reltol=1.0f-3)
381+
@trixi_test_nowarn trixi_include_changeprecision(Float32, @__MODULE__,
382+
joinpath(examples_dir(),
383+
"fluid",
384+
"pipe_flow_2d.jl"),
385+
extra_callback=steady_state_reached,
386+
dtmax=2.0f-3,
387+
tspan=(0.0f0, 1.5f0),
388+
parallelization_backend=Main.parallelization_backend,
389+
viscosity_boundary=nothing)
390+
391+
# Make sure that the simulation is terminated after a reasonable amount of time
392+
@test 0.1 < sol.t[end] < 1.0
393+
@test sol.retcode == ReturnCode.Terminated
394+
end
357395
end
358396

359397
@testset verbose=true "Solid" begin

‎test/general/custom_quantities.jl‎

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,115 @@
1+
@testset verbose=true "Custom Quantities" begin
2+
particle_spacing = 0.1
3+
coordinates = [0.0 0.1 0.2; 0.0 0.0 0.1]
4+
velocities = [1.0 2.0 0.5; 0.5 1.0 1.5]
5+
densities = [1000.0, 1000.0, 1000.0]
6+
pressures = [101325.0, 101330.0, 101320.0]
7+
8+
initial_condition = InitialCondition(; coordinates, velocity=velocities,
9+
density=densities,
10+
pressure=pressures, particle_spacing)
11+
12+
smoothing_kernel = SchoenbergCubicSplineKernel{2}()
13+
smoothing_length = 1.2 * particle_spacing
14+
15+
fluid_system = EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel,
16+
smoothing_length, 1.0)
17+
fluid_system.cache.density .= initial_condition.density
18+
19+
boundary_model = BoundaryModelDummyParticles(initial_condition.density,
20+
initial_condition.mass,
21+
AdamiPressureExtrapolation(),
22+
smoothing_kernel, smoothing_length)
23+
24+
boundary_system = BoundarySPHSystem(initial_condition, boundary_model)
25+
26+
semi = Semidiscretization(fluid_system, boundary_system)
27+
28+
v_ode, u_ode = semidiscretize(semi, (0, 1)).u0.x
29+
t = 0.0
30+
31+
@testset "Kinetic Energy" begin
32+
@testset "Fluid System" begin
33+
ekin = kinetic_energy(fluid_system, v_ode, u_ode, semi, t)
34+
35+
expected_ekin = sum(velocities) do velocity
36+
return 0.5 * first(fluid_system.mass) * dot(velocity, velocity)
37+
end
38+
39+
@test isapprox(expected_ekin, ekin)
40+
end
41+
42+
@testset "Boundary System" begin
43+
ekin = kinetic_energy(boundary_system, v_ode, u_ode, semi, t)
44+
@test ekin == 0
45+
end
46+
end
47+
48+
@testset "Total Mass" begin
49+
@testset "Fluid System" begin
50+
mass = total_mass(fluid_system, v_ode, u_ode, semi, t)
51+
expected_mass = first(fluid_system.mass) * nparticles(fluid_system)
52+
53+
@test isapprox(mass, expected_mass)
54+
end
55+
56+
@testset "Boundary System" begin
57+
mass = total_mass(boundary_system, v_ode, u_ode, semi, t)
58+
59+
@test isnan(mass)
60+
end
61+
end
62+
63+
@testset "Pressure Quantities" begin
64+
@testset "Max Pressure" begin
65+
max_p = max_pressure(fluid_system, v_ode, u_ode, semi, t)
66+
@test isapprox(max_p, 101330.0)
67+
68+
# Boundary system should return NaN
69+
@test isnan(max_pressure(boundary_system, v_ode, u_ode, semi, t))
70+
end
71+
72+
@testset "Min Pressure" begin
73+
min_p = min_pressure(fluid_system, v_ode, u_ode, semi, t)
74+
@test min_p ≈ 101320.0
75+
76+
# Boundary system should return NaN
77+
@test isnan(min_pressure(boundary_system, v_ode, u_ode, semi, t))
78+
end
79+
80+
@testset "Average Pressure" begin
81+
avg_p = avg_pressure(fluid_system, v_ode, u_ode, semi, t)
82+
expected_avg = (101325.0 + 101330.0 + 101320.0) / 3
83+
@test isapprox(avg_p, expected_avg)
84+
85+
# Boundary system should return NaN
86+
@test isnan(avg_pressure(boundary_system, v_ode, u_ode, semi, t))
87+
end
88+
end
89+
90+
@testset "Density Quantities" begin
91+
@testset "max_density" begin
92+
max_d = max_density(fluid_system, v_ode, u_ode, semi, t)
93+
@test isapprox(max_d, 1000.0) # All particles have same density
94+
95+
# Boundary system should return NaN
96+
@test isnan(max_density(boundary_system, v_ode, u_ode, semi, t))
97+
end
98+
99+
@testset "min_density" begin
100+
min_d = min_density(fluid_system, v_ode, u_ode, semi, t)
101+
@test isapprox(min_d, 1000.0)
102+
103+
# Boundary system should return NaN
104+
@test isnan(min_density(boundary_system, v_ode, u_ode, semi, t))
105+
end
106+
107+
@testset "avg_density" begin
108+
avg_d = avg_density(fluid_system, v_ode, u_ode, semi, t)
109+
@test isapprox(avg_d, 1000.0)
110+
111+
# Boundary system should return NaN
112+
@test isnan(avg_density(boundary_system, v_ode, u_ode, semi, t))
113+
end
114+
end
115+
end

‎test/general/general.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,3 +5,4 @@ include("semidiscretization.jl")
55
include("interpolation.jl")
66
include("buffer.jl")
77
include("util.jl")
8+
include("custom_quantities.jl")

0 commit comments

Comments
 (0)