BREAKING: SIMD collide forces

Project.toml

@@ -1,6 +1,6 @@
 name = "SimplexCellLists"
 uuid = "fbe69cd6-f244-4ce4-a90b-8979d86c9ea4"
-version = "0.2.3-dev"
+version = "0.3.0-dev"
 authors = ["Nathan Zimmerberg"]
 
 [workspace]

README.md

@@ -41,7 +41,7 @@ Neighbor lists and collide forces are configured with a `CollidePolicy`.
 The policy decides which objects and pairs participate in the neighbor lists
 (`filter_object`, `filter_pair`), how per-object parameters combine into
 per-pair parameters (`mix_params`), and which force law applies to each edge
-(`nl_edge_forces!`).
+(`nl_edge_force`).
 
 `DefaultCollidePolicy` stores a stiffness and collision layer masks per
 object, and applies a soft repulsive potential `E = k/2 * (L - d)²` when two
@@ -79,5 +79,5 @@ reused across steps. `setup_neighbors_naive!` is a reference implementation of
 
 A custom policy subtypes `CollidePolicy{ObjectParams, PairParams}` with its
 own parameter types and methods for `filter_object`, `filter_pair`, and
-`mix_params`, and can define its own `nl_edge_forces!` methods to change the
+`mix_params`, and can define its own `nl_edge_force` methods to change the
 force law.

benchmark/force-bench.jl

@@ -19,29 +19,44 @@ function count_interacting(pos, nl_list, policy)
     end
 end
 
+# Each system has copies of its positions, which the benchmark runs in turn.
+# When the same positions are run over and over, the branch predictor learns
+# which groups of edges have one in range, if there are fewer than about 6000
+# such groups, and the time is too low. So the ring is jittered into copies,
+# as its positions would change in a simulation. The isosphere systems have
+# enough groups already.
 systems = [
-    ("ring", Ring, Ring.pos),
-    ("isosphere", Isosphere, Isosphere.pos),
+    ("ring", Ring, [Ring.pos .+ 0.3f0 .* randn(Xoshiro(s), eltype(Ring.pos), length(Ring.pos)) for s in 1:32]),
+    ("isosphere", Isosphere, [Isosphere.pos]),
     # The isosphere is so regular that no edges are in range, so also move its vertices randomly.
-    ("isosphere jittered", Isosphere, Isosphere.pos .+ 3f0 .* randn(Xoshiro(42), eltype(Isosphere.pos), length(Isosphere.pos))),
+    ("isosphere jittered", Isosphere, [Isosphere.pos .+ 3f0 .* randn(Xoshiro(42), eltype(Isosphere.pos), length(Isosphere.pos))]),
 ]
 
-for (name, sys, pos) in systems
-    local nl = NeighborLists(sys.policy)
-    setup_neighbors_grid!(nl, pos, sys.inputs)
-    println("$name: $(length(pos)) beads")
+for (name, sys, copies) in systems
+    local nls = map(copies) do pos
+        nl = NeighborLists(sys.policy)
+        setup_neighbors_grid!(nl, pos, sys.inputs)
+        nl
+    end
+    local K = length(copies)
+    println("$name: $(length(first(copies))) beads", K > 1 ? ", $K copies, edge counts are means" : "")
     for f in (:PTNL, :CCNL, :LLNL)
-        local list = getfield(nl, f)
-        isempty(list) && continue
-        println("  $f: $(length(list)) edges, $(count_interacting(pos, list, sys.policy)) interacting")
+        isempty(getfield(first(nls), f)) && continue
+        local nedges = sum(nl -> length(getfield(nl, f)), nls)
+        local ninteracting = sum(((pos, nl),) -> count_interacting(pos, getfield(nl, f), sys.policy), zip(copies, nls))
+        println("  $f: $(round(Int, nedges/K)) edges, $(round(Int, ninteracting/K)) interacting")
     end
     for (FE, calc_type) in (
             (ForceEnergyFixedPoint{30, 30}, Float64),
             (ForceEnergyFixedPoint{30, 30}, Float32),
             (ForceEnergyFloat64, Float64),
         )
-        local fe = FE(length(pos))
+        local fe = FE(length(first(copies)))
+        # Positions in the calculation type.
+        local pos_calc = map(pos -> map(p -> SimplexCellLists.Vec3{calc_type}(p), pos), copies)
+        # Each sample is one call on the next copy.
+        local k = Ref(0)
         print("  collide_forces! $FE $calc_type:")
-        display(@b (zero_force_energy!($fe); collide_forces!($fe, $pos, $nl, $calc_type)) seconds=2)
+        display(@b (k[] = mod1(k[] + 1, $K)) (zero_force_energy!($fe); collide_forces!($fe, $pos_calc[_], $nls[_], $calc_type)) evals=1 seconds=2)
     end
 end

src/SimplexCellLists.jl

@@ -157,6 +157,7 @@ end
 include("collide-forces.jl")
 export collide_forces!
 public nl_edge_forces!
+public nl_edge_force
 
 # Precompile the collide forces for the default policy
 let NL = NeighborLists{DefaultCollidePolicy, DefaultPairParams}