Merge branch 'main' into ef/max-points-per-cell

.github/workflows/FormatCheck.yml

@@ -25,7 +25,7 @@ jobs:
         #
         # julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version = "0.13.0"))'
         run: |
-          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))'
+          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.62"))'
           julia  -e 'using JuliaFormatter; format(".")'
       - name: Format check
         run: |

benchmarks/n_body.jl

@@ -12,8 +12,15 @@ This is a more realistic benchmark for particle-based simulations than
 However, due to the higher computational cost, differences between neighborhood search
 implementations are less pronounced.
 """
-function benchmark_n_body(neighborhood_search, coordinates; parallel = true)
-    mass = 1e10 * (rand(size(coordinates, 2)) .+ 1)
+function benchmark_n_body(neighborhood_search, coordinates_; parallel = true)
+    # Passing a different backend like `CUDA.CUDABackend`
+    # allows us to change the type of the array to run the benchmark on the GPU.
+    # Passing `parallel = true` or `parallel = false` will not change anything here.
+    coordinates = PointNeighbors.Adapt.adapt(parallel, coordinates_)
+    nhs = PointNeighbors.Adapt.adapt(parallel, neighborhood_search)
+
+    # This preserves the data type of `coordinates`, which makes it work for GPU types
+    mass = 1e10 * (rand!(similar(coordinates, size(coordinates, 2))) .+ 1)
     G = 6.6743e-11
 
     dv = similar(coordinates)
@@ -36,6 +43,5 @@ function benchmark_n_body(neighborhood_search, coordinates; parallel = true)
         return dv
     end
 
-    return @belapsed $compute_acceleration!($dv, $coordinates, $mass, $G,
-                                            $neighborhood_search, $parallel)
+    return @belapsed $compute_acceleration!($dv, $coordinates, $mass, $G, $nhs, $parallel)
 end

benchmarks/plot.jl

@@ -62,7 +62,7 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations;
         neighborhood_searches = [
             TrivialNeighborhoodSearch{NDIMS}(; search_radius, eachpoint = 1:n_particles),
             GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles),
-            PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles),
+            PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles)
         ]
 
         for i in eachindex(neighborhood_searches)

benchmarks/smoothed_particle_hydrodynamics.jl

@@ -34,16 +34,80 @@ function benchmark_wcsph(neighborhood_search, coordinates; parallel = true)
                                                smoothing_length, viscosity = viscosity,
                                                density_diffusion = density_diffusion)
 
-    v = vcat(fluid.velocity, fluid.density')
-    u = copy(fluid.coordinates)
+    # Note that we cannot just disable parallelism in TrixiParticles.
+    # But passing a different backend like `CUDA.CUDABackend`
+    # allows us to change the type of the array to run the benchmark on the GPU.
+    if parallel isa Bool
+        system = fluid_system
+        nhs = neighborhood_search
+    else
+        system = PointNeighbors.Adapt.adapt(parallel, fluid_system)
+        nhs = PointNeighbors.Adapt.adapt(parallel, neighborhood_search)
+    end
+
+    v = PointNeighbors.Adapt.adapt(parallel, vcat(fluid.velocity, fluid.density'))
+    u = PointNeighbors.Adapt.adapt(parallel, coordinates)
     dv = zero(v)
 
     # Initialize the system
-    TrixiParticles.initialize!(fluid_system, neighborhood_search)
-    TrixiParticles.compute_pressure!(fluid_system, v)
+    TrixiParticles.initialize!(system, nhs)
+    TrixiParticles.compute_pressure!(system, v)
 
-    return @belapsed TrixiParticles.interact!($dv, $v, $u, $v, $u, $neighborhood_search,
-                                              $fluid_system, $fluid_system)
+    return @belapsed TrixiParticles.interact!($dv, $v, $u, $v, $u, $nhs, $system, $system)
+end
+
+"""
+    benchmark_wcsph_fp32(neighborhood_search, coordinates; parallel = true)
+
+Like [`benchmark_wcsph`](@ref), but using single precision floating point numbers.
+"""
+function benchmark_wcsph_fp32(neighborhood_search, coordinates_; parallel = true)
+    coordinates = convert(Matrix{Float32}, coordinates_)
+    density = 1000.0f0
+    fluid = InitialCondition(; coordinates, density, mass = 0.1f0)
+
+    # Compact support == smoothing length for the Wendland kernel
+    smoothing_length = convert(Float32, PointNeighbors.search_radius(neighborhood_search))
+    if ndims(neighborhood_search) == 1
+        smoothing_kernel = SchoenbergCubicSplineKernel{1}()
+    else
+        smoothing_kernel = WendlandC2Kernel{ndims(neighborhood_search)}()
+    end
+
+    sound_speed = 10.0f0
+    state_equation = StateEquationCole(; sound_speed, reference_density = density,
+                                       exponent = 1)
+
+    fluid_density_calculator = ContinuityDensity()
+    viscosity = ArtificialViscosityMonaghan(alpha = 0.02f0, beta = 0.0f0)
+    density_diffusion = DensityDiffusionMolteniColagrossi(delta = 0.1f0)
+
+    fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator,
+                                               state_equation, smoothing_kernel,
+                                               smoothing_length, viscosity = viscosity,
+                                               acceleration = (0.0f0, 0.0f0, 0.0f0),
+                                               density_diffusion = density_diffusion)
+
+    # Note that we cannot just disable parallelism in TrixiParticles.
+    # But passing a different backend like `CUDA.CUDABackend`
+    # allows us to change the type of the array to run the benchmark on the GPU.
+    if parallel isa Bool
+        system = fluid_system
+        nhs = neighborhood_search
+    else
+        system = PointNeighbors.Adapt.adapt(parallel, fluid_system)
+        nhs = PointNeighbors.Adapt.adapt(parallel, neighborhood_search)
+    end
+
+    v = PointNeighbors.Adapt.adapt(parallel, vcat(fluid.velocity, fluid.density'))
+    u = PointNeighbors.Adapt.adapt(parallel, coordinates)
+    dv = zero(v)
+
+    # Initialize the system
+    TrixiParticles.initialize!(system, nhs)
+    TrixiParticles.compute_pressure!(system, v)
+
+    return @belapsed TrixiParticles.interact!($dv, $v, $u, $v, $u, $nhs, $system, $system)
 end
 
 """

docs/make.jl

@@ -52,7 +52,7 @@ makedocs(modules = [PointNeighbors],
              "Home" => "index.md",
              "API reference" => "reference.md",
              "Authors" => "authors.md",
-             "License" => "license.md",
+             "License" => "license.md"
          ])
 
 deploydocs(;

test/neighborhood_search.jl

@@ -10,7 +10,7 @@
         names = [
             "Simple Example 2D",
             "Box Not Multiple of Search Radius 2D",
-            "Simple Example 3D",
+            "Simple Example 3D"
         ]
 
         coordinates = [
@@ -20,15 +20,15 @@
              -0.12 -0.05 -0.09 0.15 0.42],
             [-0.08 0.0 0.18 0.1 -0.08
              -0.12 -0.05 -0.09 0.15 0.39
-             0.14 0.34 0.12 0.06 0.13],
+             0.14 0.34 0.12 0.06 0.13]
         ]
 
         periodic_boxes = [
             PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.2, 0.4]),
             # The `GridNeighborhoodSearch` is forced to round up the cell sizes in this test
             # to avoid split cells.
             PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.205, 0.43]),
-            PeriodicBox(min_corner = [-0.1, -0.2, 0.05], max_corner = [0.2, 0.4, 0.35]),
+            PeriodicBox(min_corner = [-0.1, -0.2, 0.05], max_corner = [0.2, 0.4, 0.35])
         ]
 
         @testset verbose=true "$(names[i])" for i in eachindex(names)
@@ -58,15 +58,15 @@
                                                                            search_radius,
                                                                            backend = Vector{Vector{Int32}})),
                 PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points,
-                                                     periodic_box = periodic_boxes[i]),
+                                                     periodic_box = periodic_boxes[i])
             ]
 
             names = [
                 "`TrivialNeighborhoodSearch`",
                 "`GridNeighborhoodSearch`",
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors`",
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `Vector{Vector}`",
-                "`PrecomputedNeighborhoodSearch`",
+                "`PrecomputedNeighborhoodSearch`"
             ]
 
             # Also test copied templates
@@ -80,7 +80,7 @@
                                               cell_list = FullGridCellList(min_corner = periodic_boxes[i].min_corner,
                                                                            max_corner = periodic_boxes[i].max_corner,
                                                                            backend = Vector{Vector{Int32}})),
-                PrecomputedNeighborhoodSearch{NDIMS}(periodic_box = periodic_boxes[i]),
+                PrecomputedNeighborhoodSearch{NDIMS}(periodic_box = periodic_boxes[i])
             ]
             copied_nhs = copy_neighborhood_search.(template_nhs, search_radius, n_points)
             append!(neighborhood_searches, copied_nhs)
@@ -119,7 +119,7 @@
             (10, 11),
             (100, 90),
             (9, 10, 7),
-            (39, 40, 41),
+            (39, 40, 41)
         ]
 
         seeds = [1, 2]
@@ -175,7 +175,7 @@
                                                                            max_corner,
                                                                            search_radius,
                                                                            backend = Vector{Vector{Int}})),
-                PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points),
+                PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points)
             ]
 
             names = [
@@ -184,7 +184,7 @@
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors` and `ParallelUpdate`",
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors` and `SemiParallelUpdate`",
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `Vector{Vector}`",
-                "`PrecomputedNeighborhoodSearch`",
+                "`PrecomputedNeighborhoodSearch`"
             ]
 
             # Also test copied templates
@@ -199,7 +199,7 @@
                 GridNeighborhoodSearch{NDIMS}(cell_list = FullGridCellList(; min_corner,
                                                                            max_corner,
                                                                            backend = Vector{Vector{Int32}})),
-                PrecomputedNeighborhoodSearch{NDIMS}(),
+                PrecomputedNeighborhoodSearch{NDIMS}()
             ]
             copied_nhs = copy_neighborhood_search.(template_nhs, search_radius, n_points)
             append!(neighborhood_searches, copied_nhs)

test/nhs_grid.jl

@@ -200,7 +200,7 @@
         names = [
             "Simple Example 2D",
             "Box Not Multiple of Search Radius 2D",
-            "Simple Example 3D",
+            "Simple Example 3D"
         ]
 
         coordinates = [
@@ -210,15 +210,15 @@
              -0.12 -0.05 -0.09 0.15 0.42],
             [-0.08 0.0 0.18 0.1 -0.08
              -0.12 -0.05 -0.09 0.15 0.39
-             0.14 0.34 0.12 0.06 0.13],
+             0.14 0.34 0.12 0.06 0.13]
         ]
 
         periodic_boxes = [
             PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.2, 0.4]),
             # The `GridNeighborhoodSearch` is forced to round up the cell sizes in this test
             # to avoid split cells.
             PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.205, 0.43]),
-            PeriodicBox(min_corner = [-0.1, -0.2, 0.05], max_corner = [0.2, 0.4, 0.35]),
+            PeriodicBox(min_corner = [-0.1, -0.2, 0.05], max_corner = [0.2, 0.4, 0.35])
         ]
 
         @testset verbose=true "$(names[i])" for i in eachindex(names)
@@ -257,11 +257,11 @@
             nhs = GridNeighborhoodSearch{2}(search_radius = 1.0, n_points = size(coords, 2),
                                             periodic_box = PeriodicBox(min_corner = [
                                                                            -1.5,
-                                                                           0.0,
+                                                                           0.0
                                                                        ],
                                                                        max_corner = [
                                                                            2.5,
-                                                                           3.0,
+                                                                           3.0
                                                                        ]))
 
             initialize_grid!(nhs, coords)