Difference from brute force

[PharmCat]

I have data from here AhmedSalih3d
repo .

And I found different results with brute force test.

path    = dirname(@__FILE__)
DF_FLUID = CSV.read(joinpath(path, "./FluidPoints_Dp0.02.csv"), DataFrame)
DF_BOUND = CSV.read(joinpath(path, "./BoundaryPoints_Dp0.02.csv"), DataFrame)
P1F = DF_FLUID[!,"Points:0"]
P2F = DF_FLUID[!,"Points:1"]
P3F = DF_FLUID[!,"Points:2"]
P1B = DF_BOUND[!,"Points:0"]
P2B = DF_BOUND[!,"Points:1"]
P3B = DF_BOUND[!,"Points:2"]
points = Vector{SVector{3,Float64}}()
for i = 1:length(P1F)
    push!(points,SVector(P1F[i],P3F[i],P2F[i]))
end
for i = 1:length(P1B)
    push!(points,SVector(P1B[i],P3B[i],P2B[i]))
end

H = 0.06788225099390856
system  = InPlaceNeighborList(x=points, cutoff = H, parallel=true)
update!(system, points)
list = neighborlist!(system)

101128 pairs

system  = InPlaceNeighborList(x=points, cutoff = H, parallel=false)
update!(system, points)
list = neighborlist!(system)

101128 pairs

 @inline function distance_condition(p1, p2, r) 
    sum(@. (p1 - p2)^2) ≤ r^2
 end

function brute_force(p, r)
    ps = Vector{Tuple{Int, Int}}()
    n = length(p)
    for i in 1:(n-1)
        for j in (i+1):n
            if @inbounds distance_condition(p[i], p[j], r)
                    push!(ps, (i, j))
            end
        end
    end
    return ps
end

brute_force(points, H)

98585 pairs

But if change a little bit:

system  = InPlaceNeighborList(x=points, cutoff = H - 0.0000001, parallel=false)
update!(system, points)
list = neighborlist!(system)

What can I do to get more reproducible results?

Also:

brute_force(points, H + sqrt(eps()))

98585 pairs

bf = brute_force(points, H + 0.00422)

98585 pairs

bf = brute_force(points, H + 0.00423)

118425 pairs

[lmiq]

Thanks for reporting. I can reproduce the issue. There are some repeated pairs in the list coming out from neighbourlist. I'll check that asap.

[lmiq]

Thanks for reporting. This is fixed in the upcoming version 0.8.24. It was an issue associate to the fake periodic boundaries we have to add when no periodic boundary conditions are used.

I have tested the updated version with all data sets of that repository. For the records, this is the test file I'm running, after cloning the https://github.com/AhmedSalih3d/SPHExample repository. I'm running these tests from the input directory, where the files are available:

import Pkg; Pkg.activate(".")
using CSV
using CellListMap
using StaticArrays
using DataFrames
using LinearAlgebra: norm
using Test

function vprop(p, r)
    box = Box(limits(p), r)
    cl = CellList(p, box)
    vprop = map_pairwise(
        (x,y,i,j,d2,u) -> u += sqrt(d2) - r,
        0.0, box, cl
    )
    return vprop
end

function brute_force(p, r)
    ps = Vector{Tuple{Int, Int}}()
    n = length(p)
    vprop = 0.0
    for i in 1:(n-1)
        for j in (i+1):n
            d = norm(p[j] - p[i])
            if d <= r
                push!(ps, (i, j))
                vprop += d - r
            end
        end
    end
    return ps, vprop
end
    

function test(fluid_point_file, boundary_points_file)

    path    = @__DIR__
    DF_FLUID = CSV.read(joinpath(path, fluid_point_file), DataFrame)
    DF_BOUND = CSV.read(joinpath(path, boundary_points_file), DataFrame)
    
    
    P1F = DF_FLUID[!,"Points:0"]
    P2F = DF_FLUID[!,"Points:1"]
    P3F = DF_FLUID[!,"Points:2"]
    P1B = DF_BOUND[!,"Points:0"]
    P2B = DF_BOUND[!,"Points:1"]
    P3B = DF_BOUND[!,"Points:2"]
    points = Vector{SVector{3,Float64}}()
    for i = 1:length(P1F)
        push!(points,SVector(P1F[i],P3F[i],P2F[i]))
    end
    for i = 1:length(P1B)
        push!(points,SVector(P1B[i],P3B[i],P2B[i]))
    end
    
    H = 0.06788225099390856
    system  = InPlaceNeighborList(x=points, cutoff = H, parallel=true)
    update!(system, points)
    list = neighborlist!(system)
    
    list_brute_force, vprop_brute_force = brute_force(points, H)
    
    println("length(list) = ", length(list))
    println("length(list_brute_force) = ", length(list_brute_force))
    
    vprop_celllists = vprop(points, H)
    
    println("vprop by CellListMap = ", vprop_celllists)
    println("vprop by brute force= ", vprop_brute_force)
    
    @test length(list) == length(list_brute_force)
    @test vprop_celllists ≈ vprop_brute_force

    return nothing
end 
    
function run()
    files = [
        ("FluidPoints_Dp0.02.csv", "BoundaryPoints_Dp0.02.csv"),
        ("FluidPoints_Dp0.04.csv", "BoundaryPoints_Dp0.04.csv"),
        ("FluidPoints_Dp0.005.csv", "BoundaryPoints_Dp0.005.csv"),
        ("3D_DamBreak_Fluid_3LAYERS.csv", "3D_DamBreak_Boundary_3LAYERS.csv"),
        ("3D_DamBreak_Fluid.csv", "3D_DamBreak_Boundary.csv"),
    ]
    for (fluid_file, boundary_file) in files
       test(fluid_file, boundary_file)
    end
end

run()

The output is:

julia> run()
length(list) = 98585
length(list_brute_force) = 98585
vprop by CellListMap = -2143.128315646188
vprop by brute force= -2143.1283156461636
length(list) = 6635
length(list_brute_force) = 6635
vprop by CellListMap = -132.06224207331414
vprop by brute force= -132.06224207331402
length(list) = 22396347
length(list_brute_force) = 22396347
vprop by CellListMap = -514111.7302766719
vprop by brute force= -514111.7302767918
length(list) = 149377884
length(list_brute_force) = 149377884
vprop by CellListMap = -2.785453960660278e6
vprop by brute force= -2.785453960687544e6
length(list) = 123979178
length(list_brute_force) = 123979178
vprop by CellListMap = -2.2328636134960954e6
vprop by brute force= -2.2328636135205613e6

[PharmCat]

Thank you very much!

[AhmedSalih3d]

Hello guys!

Awesome you found and fixed it. Just wanted to comment it is so cool for me at a personal level to see how something I made, is being used by someone else to find a minor bug in such a cool package like CellListMap!

Made my day :)

[lmiq]

FWIW, I don't consider that a minor bug...

The issue is that all the implementation is focused on optimizing the performance for periodic systems, and the non-periodic support is an ugly hack. But it must be correct.

The coordinates of your examples are useful because they have a lot of structure, creating unusual positions when it compares to disordered systems.

[AhmedSalih3d]

FWIW, I don't consider that a minor bug...

The issue is that all the implementation is focused on optimizing the performance for periodic systems, and the non-periodic support is an ugly hack. But it must be correct.

The coordinates of your examples are useful because they have a lot of structure, creating unusual positions when it compares to disordered systems.

I apologize for poor wording, I was just very happy - I've been using CellListMap so much that I forgot it was actually made for molecular simulation and not SPH hehe 😊