Closed PharmCat closed 7 months ago
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.
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
Thank you very much!
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 :)
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.
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 😊
I have data from here AhmedSalih3d repo .
And I found different results with brute force test.
101128 pairs
101128 pairs
98585 pairs
But if change a little bit:
What can I do to get more reproducible results?
Also:
98585 pairs
98585 pairs
118425 pairs