@@ -2,19 +2,38 @@ using PointNeighbors
22using TrixiParticles
33using BenchmarkTools
44
5+ # Create a dummy semidiscretization type to be able to use a specific neighborhood search
6+ struct DummySemidiscretization{N, P}
7+ neighborhood_search :: N
8+ parallelization_backend :: P
9+ end
10+
11+ @inline function PointNeighbors. parallel_foreach (f, iterator, semi:: DummySemidiscretization )
12+ PointNeighbors. parallel_foreach (f, iterator, semi. parallelization_backend)
13+ end
14+
15+ @inline function TrixiParticles. get_neighborhood_search (_, _, semi:: DummySemidiscretization )
16+ return semi. neighborhood_search
17+ end
18+
19+ @inline function TrixiParticles. get_neighborhood_search (_, semi:: DummySemidiscretization )
20+ return semi. neighborhood_search
21+ end
22+
523"""
624 benchmark_wcsph(neighborhood_search, coordinates; parallel = true)
725
826A benchmark of the right-hand side of a full real-life Weakly Compressible
927Smoothed Particle Hydrodynamics (WCSPH) simulation with TrixiParticles.jl.
1028This method is used to simulate an incompressible fluid.
1129"""
12- function benchmark_wcsph (neighborhood_search, coordinates; parallel = true )
30+ function benchmark_wcsph (neighborhood_search, coordinates;
31+ parallelization_backend = default_backend (coordinates))
1332 density = 1000.0
1433 fluid = InitialCondition (; coordinates, density, mass = 0.1 )
1534
16- # Compact support == smoothing length for the Wendland kernel
17- smoothing_length = PointNeighbors. search_radius (neighborhood_search)
35+ # Compact support == 2 * smoothing length for these kernels
36+ smoothing_length = PointNeighbors. search_radius (neighborhood_search) / 2
1837 if ndims (neighborhood_search) == 1
1938 smoothing_kernel = SchoenbergCubicSplineKernel {1} ()
2039 else
@@ -34,40 +53,36 @@ function benchmark_wcsph(neighborhood_search, coordinates; parallel = true)
3453 smoothing_length, viscosity = viscosity,
3554 density_diffusion = density_diffusion)
3655
37- # Note that we cannot just disable parallelism in TrixiParticles.
38- # But passing a different backend like `CUDA.CUDABackend`
39- # allows us to change the type of the array to run the benchmark on the GPU.
40- if parallel isa Bool
41- system = fluid_system
42- nhs = neighborhood_search
43- else
44- system = PointNeighbors. Adapt. adapt (parallel, fluid_system)
45- nhs = PointNeighbors. Adapt. adapt (parallel, neighborhood_search)
46- end
56+ system = PointNeighbors. Adapt. adapt (parallelization_backend, fluid_system)
57+ nhs = PointNeighbors. Adapt. adapt (parallelization_backend, neighborhood_search)
58+ semi = DummySemidiscretization (nhs, parallelization_backend)
4759
48- v = PointNeighbors. Adapt. adapt (parallel, vcat (fluid. velocity, fluid. density' ))
49- u = PointNeighbors. Adapt. adapt (parallel, coordinates)
60+ v = PointNeighbors. Adapt. adapt (parallelization_backend,
61+ vcat (fluid. velocity, fluid. density' ))
62+ u = PointNeighbors. Adapt. adapt (parallelization_backend, coordinates)
5063 dv = zero (v)
5164
5265 # Initialize the system
53- TrixiParticles. initialize! (system, nhs )
54- TrixiParticles. compute_pressure! (system, v)
66+ TrixiParticles. initialize! (system, semi )
67+ TrixiParticles. compute_pressure! (system, v, semi )
5568
56- return @belapsed TrixiParticles. interact! ($ dv, $ v, $ u, $ v, $ u, $ nhs , $ system, $ system )
69+ return @belapsed TrixiParticles. interact! ($ dv, $ v, $ u, $ v, $ u, $ system , $ system, $ semi )
5770end
5871
5972"""
6073 benchmark_wcsph_fp32(neighborhood_search, coordinates; parallel = true)
6174
6275Like [`benchmark_wcsph`](@ref), but using single precision floating point numbers.
6376"""
64- function benchmark_wcsph_fp32 (neighborhood_search, coordinates_; parallel = true )
77+ function benchmark_wcsph_fp32 (neighborhood_search, coordinates_;
78+ parallelization_backend = default_backend (coordinates_))
6579 coordinates = convert (Matrix{Float32}, coordinates_)
6680 density = 1000.0f0
6781 fluid = InitialCondition (; coordinates, density, mass = 0.1f0 )
6882
69- # Compact support == smoothing length for the Wendland kernel
70- smoothing_length = convert (Float32, PointNeighbors. search_radius (neighborhood_search))
83+ # Compact support == 2 * smoothing length for these kernels
84+ smoothing_length = convert (Float32,
85+ PointNeighbors. search_radius (neighborhood_search) / 2 )
7186 if ndims (neighborhood_search) == 1
7287 smoothing_kernel = SchoenbergCubicSplineKernel {1} ()
7388 else
@@ -85,29 +100,24 @@ function benchmark_wcsph_fp32(neighborhood_search, coordinates_; parallel = true
85100 fluid_system = WeaklyCompressibleSPHSystem (fluid, fluid_density_calculator,
86101 state_equation, smoothing_kernel,
87102 smoothing_length, viscosity = viscosity,
88- acceleration = (0.0f0 , 0.0f0 , 0.0f0 ),
103+ acceleration = ntuple (_ -> 0.0f0 ,
104+ Val (ndims (neighborhood_search))),
89105 density_diffusion = density_diffusion)
90106
91- # Note that we cannot just disable parallelism in TrixiParticles.
92- # But passing a different backend like `CUDA.CUDABackend`
93- # allows us to change the type of the array to run the benchmark on the GPU.
94- if parallel isa Bool
95- system = fluid_system
96- nhs = neighborhood_search
97- else
98- system = PointNeighbors. Adapt. adapt (parallel, fluid_system)
99- nhs = PointNeighbors. Adapt. adapt (parallel, neighborhood_search)
100- end
107+ system = PointNeighbors. Adapt. adapt (parallelization_backend, fluid_system)
108+ nhs = PointNeighbors. Adapt. adapt (parallelization_backend, neighborhood_search)
109+ semi = DummySemidiscretization (nhs, parallelization_backend)
101110
102- v = PointNeighbors. Adapt. adapt (parallel, vcat (fluid. velocity, fluid. density' ))
103- u = PointNeighbors. Adapt. adapt (parallel, coordinates)
111+ v = PointNeighbors. Adapt. adapt (parallelization_backend,
112+ vcat (fluid. velocity, fluid. density' ))
113+ u = PointNeighbors. Adapt. adapt (parallelization_backend, coordinates)
104114 dv = zero (v)
105115
106116 # Initialize the system
107- TrixiParticles. initialize! (system, nhs )
108- TrixiParticles. compute_pressure! (system, v)
117+ TrixiParticles. initialize! (system, semi )
118+ TrixiParticles. compute_pressure! (system, v, semi )
109119
110- return @belapsed TrixiParticles. interact! ($ dv, $ v, $ u, $ v, $ u, $ nhs , $ system, $ system )
120+ return @belapsed TrixiParticles. interact! ($ dv, $ v, $ u, $ v, $ u, $ system , $ system, $ semi )
111121end
112122
113123"""
@@ -117,12 +127,13 @@ A benchmark of the right-hand side of a full real-life Total Lagrangian
117127Smoothed Particle Hydrodynamics (TLSPH) simulation with TrixiParticles.jl.
118128This method is used to simulate an elastic structure.
119129"""
120- function benchmark_tlsph (neighborhood_search, coordinates; parallel = true )
130+ function benchmark_tlsph (neighborhood_search, coordinates;
131+ parallelization_backend = default_backend (coordinates))
121132 material = (density = 1000.0 , E = 1.4e6 , nu = 0.4 )
122133 solid = InitialCondition (; coordinates, density = material. density, mass = 0.1 )
123134
124- # Compact support == smoothing length for the Wendland kernel
125- smoothing_length = PointNeighbors. search_radius (neighborhood_search)
135+ # Compact support == 2 * smoothing length for these kernels
136+ smoothing_length = PointNeighbors. search_radius (neighborhood_search) / 2
126137 if ndims (neighborhood_search) == 1
127138 smoothing_kernel = SchoenbergCubicSplineKernel {1} ()
128139 else
@@ -131,14 +142,15 @@ function benchmark_tlsph(neighborhood_search, coordinates; parallel = true)
131142
132143 solid_system = TotalLagrangianSPHSystem (solid, smoothing_kernel, smoothing_length,
133144 material. E, material. nu)
145+ semi = DummySemidiscretization (neighborhood_search, parallelization_backend)
134146
135147 v = copy (solid. velocity)
136148 u = copy (solid. coordinates)
137149 dv = zero (v)
138150
139151 # Initialize the system
140- TrixiParticles. initialize! (solid_system, neighborhood_search )
152+ TrixiParticles. initialize! (solid_system, semi )
141153
142- return @belapsed TrixiParticles. interact! ($ dv, $ v, $ u, $ v, $ u, $ neighborhood_search,
143- $ solid_system, $ solid_system)
154+ return @belapsed TrixiParticles. interact! ($ dv, $ v, $ u, $ v, $ u,
155+ $ solid_system, $ solid_system, $ semi )
144156end
0 commit comments