Bad scaling when 2 sets of particles are used.

[lmiq]

See: trixi-framework/PointNeighbors.jl#8 (comment)

[lmiq]

The issue, as explained in trixi-framework/PointNeighbors.jl#8 (comment) is that when using 2 sets of particles we are not using the projection-based approach to avoid computing distances for particles that are not within the cutoff.

As a very early test, here is an implementation in which we just duplicate the particles and construct the cell lists for all duplicated set, and use the single-set approach. This is far from ideal, since we are still looping over much more particles than needed, but nevertheless it scales much better already:

julia> run(;n=10^5)
CellListMap (t1): 0.097850265 s
CellListMap different (t2/t1): 1.755102543667102
CellListMap single_list (t5/t1): 2.697372194137645
Test Passed

julia> run(;n=5*10^5)
CellListMap (t1): 0.45006516 s
CellListMap different (t2/t1): 2.006997344562285
CellListMap single_list (t5/t1): 2.740345904579683
Test Passed

julia> run(;n=10^6)
CellListMap (t1): 0.90805197 s
CellListMap different (t2/t1): 2.940807859268231
CellListMap single_list (t5/t1): 2.605681498604094
Test Passed

julia> run(;n=5*10^6)
CellListMap (t1): 4.56837243 s
CellListMap different (t2/t1): 5.094499038906073
CellListMap single_list (t5/t1): 2.7785865216772616
Test Passed

Note that here we are using a triclinic cell, thus ideally we should expect that all runs were comparable to t1 (no 2x overhead implied - we are not using symmetry). These benchmarks to do not include initialization of the ParticleSystems, where the overhead of creating cell lists for a greater number of particles would appear. These are the initialization times for different sizes of systems:

julia> x, box = CellListMap.xatomic(10^5);

julia> @b ParticleSystem(positions=x, cutoff=box.cutoff, output=similar(x), unitcell=box.input_unit_cell.matrix)
16.301 ms (47290 allocs: 92.691 MiB)

julia> x, box = CellListMap.xatomic(5*10^5);

julia> @b ParticleSystem(positions=x, cutoff=box.cutoff, output=similar(x), unitcell=box.input_unit_cell.matrix)
107.122 ms (139089 allocs: 403.936 MiB, 28.82% gc time)

julia> x, box = CellListMap.xatomic(10^6);

julia> @b ParticleSystem(positions=x, cutoff=box.cutoff, output=similar(x), unitcell=box.input_unit_cell.matrix)
209.146 ms (237622 allocs: 789.580 MiB, 18.13% gc time, without a warmup)

julia> x, box = CellListMap.xatomic(5*10^6);

julia> @b ParticleSystem(positions=x, cutoff=box.cutoff, output=similar(x), unitcell=box.input_unit_cell.matrix)
1.685 s (1036615 allocs: 3.685 GiB, 31.67% gc time, without a warmup)

Which, we see, are not negligible. The fraction of time spent on initialization relative to the time spent on mapping, in these examples, are, thus:

10^5       :    0.16
5*10^5   :    0.23
10^6       :    0.23
5*10^6  :     0.36

Remembering that this initialization time is largely dependent on the size of the smallest set, reason why it is avoided now:

julia> x, box = CellListMap.xatomic(5*10^6);

julia> xypositions = vcat(x,y);

julia> @b ParticleSystem(positions=xypositions, cutoff=box.cutoff, output=similar(x), unitcell=box.input_unit_cell.matrix)
2.381 s (1322650 allocs: 4.989 GiB, 35.57% gc time, without a warmup)

julia> y = x[1:10^4];

julia> @b ParticleSystem(xpositions=x, ypositions=y, cutoff=box.cutoff, output=similar(x), unitcell=box.input_unit_cell.matrix)
981.499 ms (46537 allocs: 2.475 GiB, 29.04% gc time, without a warmup)

Code:

Details

using CellListMap, PointNeighbors, StaticArrays
using Chairmarks
using Test

function test_celllistmap(_, sys)
    sys.output .= 0
    map_pairwise!(sys) do x, y, i, j, d2, f
        d = sqrt(d2)
        fij = (y - x) / d^3
        for dim in axes(f, 1)
            @inbounds f[dim, i] += fij[dim]
            @inbounds f[dim, j] -= fij[dim]
        end

        return f
    end
    return sys.output
end

function test_celllistmap_single_list(_, sys)
    sys.output .= 0
    map_pairwise!(sys) do x, y, i, j, d2, f
        d2 < eps() && return f
        n = size(f,2)
        if (i <= n && j > n) ||
           (j <= n && i > n)
            d = sqrt(d2)
            if i < j
                fij = (y - x) / d^3
            else
                fij = (x - y) / d^3
            end
            for dim in axes(f, 1)
                f[dim, min(i,j)] += fij[dim]
            end
        end
        return f
    end
    return sys.output
end

function test_celllistmap_different(_, sys)
    sys.output .= 0
    map_pairwise!(sys) do x, y, i, j, d2, f
        d2 < eps() && return f
        d = sqrt(d2)
        fij = (y - x) / d^3
        for dim in axes(f, 1)
            @inbounds f[dim, i] += fij[dim]
        end
        return f
    end
    return sys.output
end

function test_pointneighbors(x, f, nhs)
    f .= 0
    foreach_point_neighbor(x, x, nhs, parallel=true) do i, j, pos_diff, distance
        distance < sqrt(eps()) && return
        dv = -pos_diff / distance^3
        for dim in axes(f, 1)
            @inbounds f[dim, i] += dv[dim]
        end
    end
    return f
end

function run(;n=10^5)
    x, box = CellListMap.xatomic(n);
    uc = box.input_unit_cell.matrix

    xmat = Matrix(reinterpret(reshape, Float64, x));
    f = similar(xmat)
    f2 = similar(xmat)
    f3 = similar(xmat)
    f4 = similar(xmat)
    f5 = similar(xmat)

    xymat = hcat(xmat, xmat)

    nhs1 = GridNeighborhoodSearch{3}(; search_radius=box.cutoff, n_points=size(xmat,2))
    initialize!(nhs1, xmat, xmat)
    
    min_corner = minimum(xmat, dims=2) .- 5
    max_corner = maximum(xmat, dims=2) .+ 5
    nhs2 = GridNeighborhoodSearch{3}(; search_radius=box.cutoff, n_points=size(xmat,2),
                                    cell_list=FullGridCellList(search_radius=box.cutoff;
                                                                min_corner, max_corner,
                                                                max_points_per_cell=2000))
    initialize!(nhs2, xmat, xmat)
    sys = ParticleSystem(positions=x, cutoff=box.cutoff, output=f, unitcell=uc)
    sys2 = ParticleSystem(xpositions=x, ypositions=x, cutoff=box.cutoff, output=f2, unitcell=uc)
    sys3 = ParticleSystem(xpositions=xymat, cutoff=box.cutoff, output=f5, unitcell=uc)

    b1 = @b test_celllistmap($(nothing), $sys)
    println("CellListMap (t1): ", b1.time, " s")
    b2 = @b test_celllistmap_different($(nothing), $sys2)
    println("CellListMap different (t2/t1): ", b2.time / b1.time) 
    b5 = @b test_celllistmap_single_list($(nothing), $sys3)
    println("CellListMap single_list (t5/t1): ", b5.time / b1.time) 
#    b3 = @b test_pointneighbors($xmat, $f3, $nhs1)
#    println("PointNeighbors nhs1 (t3/t1): ", b3.time / b1.time)
#    b4 = @b test_pointneighbors($xmat, $f4, $nhs2)
#    println("PointNeighbors nhs2 (t4/t1): ", b4.time / b1.time)

#    @test f ≈ f2 ≈ f3 ≈ f4
    @test f ≈ f2 ≈ f5
end