Merge branch 'main' into ef/neighbor-lists-vectorofvectors

Author: efaulhaber

.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="2.1.2"))'
           julia  -e 'using JuliaFormatter; format(".")'
       - name: Format check
         run: |

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/typos@v1.21.0
+        uses: crate-ci/typos@v1.31.2

.github/workflows/ci.yml

@@ -34,58 +34,69 @@ concurrency:
 
 jobs:
   build:
-    name: Run Tests (Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }})
-    runs-on: ubuntu-latest
+    name: Run Tests (Julia ${{ matrix.version }} - ${{ matrix.os }})
+    runs-on: ${{ matrix.os }}
     strategy:
       # Don't cancel all running jobs when one job fails
       fail-fast: false
       matrix:
         version:
-          - '1.9'
-          - '1'
+          - 'min' # the earliest supported version compatible with the project
+          - '1'   # the latest stable 1.x release
         os:
           - ubuntu-latest
-        arch:
-          - x64
-        include:
-          # Also run tests on Windows and macOS-ARM, but only with the latest Julia version
-          - version: '1'
-            os: windows-latest
-            arch: x64
-          - version: '1'
-            os: macos-14
-            arch: arm64
+          - windows-latest
+          - macos-14
 
     steps:
       - name: Check out project
         uses: actions/checkout@v4
+
       - name: Set up Julia
         uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
-      - run: julia -e 'using InteractiveUtils; versioninfo(verbose=true)'
+
+      - name: Display Julia version
+        run: julia -e 'using InteractiveUtils; versioninfo(verbose=true)'
+
       - uses: julia-actions/cache@v2
+
       - name: Build package
         uses: julia-actions/julia-buildpkg@v1
-      - name: Run tests
+
+      - name: Run unit tests
         uses: julia-actions/julia-runtest@v1
         with:
           annotate: true
           # Only run coverage in one Job (Ubuntu and latest Julia version)
           coverage: ${{ matrix.os == 'ubuntu-latest' && matrix.version == '1' }}
+        env:
+          POINTNEIGHBORS_TEST: unit
+
       - name: Process coverage results
         # Only run coverage in one Job (Ubuntu and latest Julia version)
         if: matrix.os == 'ubuntu-latest' && matrix.version == '1'
         uses: julia-actions/julia-processcoverage@v1
         with:
           directories: src,test
+
       - name: Upload coverage report to Codecov
         # Only run coverage in one Job (Ubuntu and latest Julia version)
         if: matrix.os == 'ubuntu-latest' && matrix.version == '1'
-        uses: codecov/codecov-action@v4
+        uses: codecov/codecov-action@v5
         with:
           files: lcov.info
           fail_ci_if_error: true
           flags: unit
         env:
           CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
+
+      - name: Run benchmark tests
+        uses: julia-actions/julia-runtest@v1
+        with:
+          annotate: true
+          coverage: false
+        env:
+          POINTNEIGHBORS_TEST: benchmarks
+

.gitignore

@@ -9,6 +9,7 @@ coverage_report/
 *.jl.*.cov
 .vscode/
 run
+out/*
 
 .DS_Store
 

.zenodo.json

@@ -0,0 +1,34 @@
+{
+    "description": "Neighborhood search with fixed search radius in Julia",
+    "license": "MIT",
+    "title": "PointNeighbors.jl",
+    "upload_type": "software",
+    "creators": [
+        {
+            "affiliation": "Department of Mathematics and Computer Science, University of Cologne, Germany",
+            "name": "Faulhaber, Erik",
+            "orcid": "0000-0001-9788-5949"
+        },
+        {
+            "affiliation": "High-Performance Computing Center Stuttgart, University of Stuttgart, Germany",
+            "name": "Neher, Niklas",
+            "orcid": "0009-0004-2472-0923"
+        },
+        {
+            "affiliation": "Institute of Surface Science, Helmholtz-Zentrum hereon, Germany",
+            "name": "Berger, Sven",
+            "orcid": "0000-0001-6083-7038"
+        },
+        {
+            "affiliation": "High-Performance Scientific Computing, University of Augsburg, Germany",
+            "name": "Schlottke-Lakemper, Michael",
+            "orcid": "0000-0002-3195-2536"
+        },
+        {
+            "affiliation": "Department of Mathematics and Computer Science, University of Cologne, Germany",
+            "name": "Gassner, Gregor J.",
+            "orcid": "0000-0002-1752-1158"
+        }
+    ],
+    "access_right": "open"
+}

Project.toml

@@ -1,19 +1,27 @@
 name = "PointNeighbors"
 uuid = "1c4d5385-0a27-49de-8e2c-43b175c8985c"
 authors = ["Erik Faulhaber <[email protected]>"]
-version = "0.2.4-dev"
+version = "0.6.3"
 
 [deps]
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ArraysOfArrays = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018"
+Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
+GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
+Adapt = "4"
 ArraysOfArrays = "0.6"
+Atomix = "1"
+GPUArraysCore = "0.2"
+KernelAbstractions = "0.9"
 LinearAlgebra = "1"
 Polyester = "0.7.5"
 Reexport = "1"
 StaticArrays = "1"
-julia = "1.9"
+julia = "1.10"

README.md

@@ -8,5 +8,63 @@
 [![Codecov](https://codecov.io/gh/trixi-framework/PointNeighbors.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/trixi-framework/PointNeighbors.jl)
 [![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle)
 [![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/license/mit/)
+[![DOI](https://zenodo.org/badge/doi/10.5281/zenodo.12702157.svg)](https://zenodo.org/doi/10.5281/zenodo.12702157)
 
-Work in Progress!
+**PointNeighbors.jl** is a package for neighborhood search with fixed search radius in
+1D, 2D and 3D point clouds.
+
+## Features
+
+- Several implementations of neighborhood search with fixed search radius
+- Focus on fast incremental updates to be usable for particle-based simulations with
+  frequent updates
+- Designed as a "playground" to easily switch between different implementations and data
+  structures
+- Common API over all implementations
+- Extensive benchmark suite to study different implementations (work in progress)
+- GPU compatibility (work in progress)
+
+| Implementation  | Description | Features | Query | Update | GPU-compatible |
+| ------------- | ------------- | --- | :--: | :--: | :--: |
+| `GridNeighborhoodSearch` with `DictionaryCellList` | Grid-based NHS with Julia `Dict` backend | Infinite domain | Fast | Fast | ❌ |
+| `GridNeighborhoodSearch` with `FullGridCellList` | Grid-based NHS allocating all cells of the domain | Finite domain, but efficient memory layout for densely filled domain. | Faster | Fastest | ✅ |
+| `PrecomputedNeighborhoodSearch` | Precompute neighbor lists | Best for [TLSPH](https://trixi-framework.github.io/TrixiParticles.jl/stable/systems/total_lagrangian_sph/) without NHS updates. Not suitable for updates in every time step. | Fastest | Very slow | ❌ |
+
+## Benchmarks
+
+The following benchmarks were conducted on an AMD Ryzen Threadripper 3990X using 128 threads.
+
+Benchmark of a single force computation step of a Weakly Compressible SPH (WCSPH) simulation:
+![wcsph](https://github.com/trixi-framework/PointNeighbors.jl/assets/44124897/ad5c378b-9ce2-4e6f-91dc-1e0da379b91f)
+
+Benchmark of an incremental update similar to a WCSPH simulation (note the log scale):
+![update](https://github.com/trixi-framework/PointNeighbors.jl/assets/44124897/71eac5c9-6aa5-4267-bc0b-4057c89f8b12)
+
+Benchmark of a full right-hand side evaluation of a WCSPH simulation (note the log scale):
+![rhs](https://github.com/trixi-framework/PointNeighbors.jl/assets/44124897/ac328a96-1b9f-4319-a785-dce9d862fd70)
+
+
+## Packages using PointNeighbors.jl
+
+- [TrixiParticles.jl](https://github.com/trixi-framework/TrixiParticles.jl)
+- [Peridynamics.jl](https://github.com/kaipartmann/Peridynamics.jl)
+- [PeriLab.jl](https://github.com/PeriHub/PeriLab.jl)
+
+If you're using PointNeighbors.jl in your package, please feel free to open a PR adding it
+to this list.
+
+
+## Cite Us
+
+If you use PointNeighbors.jl in your own research or write a paper using results obtained
+with the help of PointNeighbors.jl, please cite it as
+```bibtex
+@misc{pointneighbors,
+  title={{P}oint{N}eighbors.jl: {N}eighborhood search with fixed search radius in {J}ulia},
+  author={Erik Faulhaber and Niklas Neher and Sven Berger and
+          Michael Schlottke-Lakemper and Gregor Gassner},
+  year={2024},
+  howpublished={\url{https://github.com/trixi-framework/PointNeighbors.jl}},
+  doi={10.5281/zenodo.12702157}
+}
+```

benchmarks/benchmarks.jl

@@ -0,0 +1,6 @@
+include("count_neighbors.jl")
+include("n_body.jl")
+include("smoothed_particle_hydrodynamics.jl")
+include("update.jl")
+
+include("plot.jl")

benchmarks/count_neighbors.jl

@@ -0,0 +1,31 @@
+using PointNeighbors
+using BenchmarkTools
+
+"""
+    benchmark_count_neighbors(neighborhood_search, coordinates; parallel = true)
+
+A very cheap and simple neighborhood search benchmark, only counting the neighbors of each
+point. For each point-neighbor pair, only an array entry is incremented.
+
+Due to the minimal computational cost, differences between neighborhood search
+implementations are highlighted. On the other hand, this is the least realistic benchmark.
+
+For a computationally heavier benchmark, see [`benchmark_n_body`](@ref).
+"""
+function benchmark_count_neighbors(neighborhood_search, coordinates;
+                                   parallelization_backend = default_backend(coordinates))
+    n_neighbors = zeros(Int, size(coordinates, 2))
+
+    function count_neighbors!(n_neighbors, coordinates, neighborhood_search,
+                              parallelization_backend)
+        n_neighbors .= 0
+
+        foreach_point_neighbor(coordinates, coordinates, neighborhood_search;
+                               parallelization_backend) do i, _, _, _
+            n_neighbors[i] += 1
+        end
+    end
+
+    return @belapsed $count_neighbors!($n_neighbors, $coordinates,
+                                       $neighborhood_search, $parallelization_backend)
+end

benchmarks/n_body.jl

@@ -0,0 +1,50 @@
+using PointNeighbors
+using BenchmarkTools
+
+"""
+    benchmark_n_body(neighborhood_search, coordinates; parallel = true)
+
+A simple neighborhood search benchmark, computing the right-hand side of an n-body
+simulation with a cutoff (corresponding to the search radius of `neighborhood_search`).
+
+This is a more realistic benchmark for particle-based simulations than
+[`benchmark_count_neighbors`](@ref).
+However, due to the higher computational cost, differences between neighborhood search
+implementations are less pronounced.
+"""
+function benchmark_n_body(neighborhood_search, coordinates_;
+                          parallelization_backend = default_backend(coordinates_))
+    # 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(parallelization_backend, coordinates_)
+    nhs = PointNeighbors.Adapt.adapt(parallelization_backend, 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)
+
+    function compute_acceleration!(dv, coordinates, mass, G, neighborhood_search,
+                                   parallelization_backend)
+        dv .= 0.0
+
+        foreach_point_neighbor(coordinates, coordinates, neighborhood_search;
+                               parallelization_backend) do i, j, pos_diff, distance
+            # Only consider particles with a distance > 0
+            distance < sqrt(eps()) && return
+
+            dv_ = -G * mass[j] * pos_diff / distance^3
+
+            for dim in axes(dv, 1)
+                @inbounds dv[dim, i] += dv_[dim]
+            end
+        end
+
+        return dv
+    end
+
+    return @belapsed $compute_acceleration!($dv, $coordinates, $mass, $G, $nhs,
+                                            $parallelization_backend)
+end