From f6e71a542ef9cb93ad51317573b4eece724db999 Mon Sep 17 00:00:00 2001 From: nhz2 Date: Mon, 10 Aug 2026 09:46:25 +0200 Subject: [PATCH 1/3] add force energy --- src/SimplexCellLists.jl | 22 +++ src/force-energy.jl | 329 ++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 + test/test-force-energy.jl | 292 +++++++++++++++++++++++++++++++++ test/test-mindistance.jl | 1 + 5 files changed, 646 insertions(+) create mode 100644 src/force-energy.jl create mode 100644 test/test-force-energy.jl diff --git a/src/SimplexCellLists.jl b/src/SimplexCellLists.jl index 2845075..ce85b80 100644 --- a/src/SimplexCellLists.jl +++ b/src/SimplexCellLists.jl @@ -46,4 +46,26 @@ let T = Int64, F = Float32 precompile(cell_line_seg_add!, (LineSegCellList{T,F}, SVector{3,F}, SVector{3,F}, T)) end +include("force-energy.jl") +export ForceEnergyT +export ForceEnergyFloat64 +export ForceEnergyFloat32 +export ForceEnergyFixedPoint +export ForceNoEnergyFixedPoint +export DebugForceEnergy +export NullForceEnergy +export add_bead_force! +export zero_bead_force! +export add_energy! +export zero_energy! +export get_nbeads +export zero_force_energy! +export get_force! +export get_bead_force +export get_energy +export get_force +export combine_force_energy! + + + end \ No newline at end of file diff --git a/src/force-energy.jl b/src/force-energy.jl new file mode 100644 index 0000000..248906f --- /dev/null +++ b/src/force-energy.jl @@ -0,0 +1,329 @@ +""" + abstract type ForceEnergy + +Accumulator interface for summing per-bead forces and total energy during a +force calculation. Concrete strategies trade off precision, thread +combinability, and whether energy is tracked at all — see +[`ForceEnergyT`](@ref), [`ForceEnergyFixedPoint`](@ref), +[`ForceNoEnergyFixedPoint`](@ref), [`DebugForceEnergy`](@ref), and +[`NullForceEnergy`](@ref). + +# Implementing a custom accumulator + +A concrete `FE <: ForceEnergy` provides a constructor `FE(nbeads)` and the +write/read methods below. + +Required, called in the hot force loop: + - [`add_bead_force!`](@ref)`(fe, i, force)` — add an `SVector{3}` to bead `i`. + Annotate `Base.@propagate_inbounds`. + - [`zero_bead_force!`](@ref)`(fe, i)` — set force of bead `i` to zero. Useful for pinning. + Annotate `Base.@propagate_inbounds`. + - [`add_energy!`](@ref)`(fe, energy)` — add a scalar to the total energy. + - [`zero_energy!`](@ref)`(fe)` — set the total energy to zero. + +Required, for setup and reading results (accumulators that discard data, like +[`NullForceEnergy`](@ref) and [`ForceNoEnergyFixedPoint`](@ref), leave +the readers for the discarded data unimplemented so reading errors instead of +returning made-up values): + - [`get_nbeads`](@ref)`(fe)` — number of beads. + - [`zero_force_energy!`](@ref)`(fe)` — reset all forces and energy to zero. + - [`get_force!`](@ref)`(fe, out)` — write per-bead `SVector{3}` forces into `out`. + - [`get_bead_force`](@ref)`(fe, i, F=Float64)` — one bead's force. + - [`get_energy`](@ref)`(fe, F=Float64)` — the accumulated energy. + +Provided generically from the above, and may be specialized for speed or +bit-exactness: [`get_force`](@ref), +[`combine_force_energy!`](@ref) (merge one accumulator into another of matching +size, e.g. to reduce per-thread partials). +""" +abstract type ForceEnergy end + +function get_force(force_energy::ForceEnergy, F=Float64) + force_out = zeros(SVector{3, F}, get_nbeads(force_energy)) + get_force!(force_energy, force_out) + force_out +end +function combine_force_energy!(fe_out::ForceEnergy, fe_in::ForceEnergy) + @argcheck get_nbeads(fe_in) == get_nbeads(fe_out) + add_energy!(fe_out, get_energy(fe_in)) + for i in 1:get_nbeads(fe_in) + @inbounds add_bead_force!(fe_out, i, get_bead_force(fe_in, i)) + end + fe_out +end + +mutable struct ForceEnergyT{T} <: ForceEnergy + const forces::Memory{SVector{3,T}} + energy::T +end +const ForceEnergyFloat64 = ForceEnergyT{Float64} +const ForceEnergyFloat32 = ForceEnergyT{Float32} +function ForceEnergyT{T}(nbeads::Integer) where T + forces = Memory{SVector{3,T}}(undef, nbeads) + fill!(forces, zero(SVector{3,T})) + ForceEnergyT{T}(forces, zero(T)) +end +function zero_force_energy!(force_energy::ForceEnergyT{T}) where T + fill!(force_energy.forces, zero(SVector{3, T})) + force_energy.energy = zero(T) + nothing +end +function get_nbeads(force_energy::ForceEnergyT) + length(force_energy.forces) +end +function add_energy!(force_energy::ForceEnergyT, energy) + force_energy.energy += energy + nothing +end +function zero_energy!(force_energy::ForceEnergyT{T}) where {T} + force_energy.energy = zero(T) + nothing +end +Base.@propagate_inbounds function add_bead_force!(force_energy::ForceEnergyT, i::Integer, force) + force_energy.forces[i] += force + nothing +end +Base.@propagate_inbounds function zero_bead_force!(force_energy::ForceEnergyT{T}, i::Integer) where {T} + force_energy.forces[i] = zero(SVector{3,T}) + nothing +end +function get_force!(force_energy::ForceEnergyT, force_out) + Base.require_one_based_indexing(force_out) + @argcheck length(force_out) == length(force_energy.forces) + @inbounds for i in eachindex(force_out) + force_out[i] = force_energy.forces[i] + end + nothing +end +Base.@propagate_inbounds function get_bead_force(force_energy::ForceEnergyT{T}, i::Integer, ::Type{F}=Float64)::SVector{3, F} where {T, F} + SVector{3, F}(force_energy.forces[i]) +end +function get_energy(force_energy::ForceEnergyT, ::Type{F}=Float64) where F + F(force_energy.energy) +end + + +@inline convert_to_fixed_point(x, fbits) = unsafe_trunc(Int64, x*(oftype(x, 2)^fbits)) +@inline convert_to_floating_point(T, x, fbits) = T(x)*(convert(T, 2)^-fbits) + +mutable struct ForceEnergyFixedPoint{F_FBITS, E_FBITS} <: ForceEnergy + const forces::Memory{Int64} + energy::Int64 +end +function ForceEnergyFixedPoint{F_FBITS, E_FBITS}(nbeads::Integer) where {F_FBITS, E_FBITS} + forces = Memory{Int64}(undef, 3*nbeads) + fill!(forces, Int64(0)) + ForceEnergyFixedPoint{F_FBITS, E_FBITS}(forces, Int64(0)) +end +function zero_force_energy!(force_energy::ForceEnergyFixedPoint) + fill!(force_energy.forces, Int64(0)) + force_energy.energy = 0 + nothing +end +function get_nbeads(force_energy::ForceEnergyFixedPoint) + length(force_energy.forces)÷3 +end +function add_energy!(force_energy::ForceEnergyFixedPoint{F_FBITS, E_FBITS}, energy) where {F_FBITS, E_FBITS} + force_energy.energy += convert_to_fixed_point(energy, E_FBITS) + nothing +end +function zero_energy!(force_energy::ForceEnergyFixedPoint) + force_energy.energy = 0 + nothing +end +Base.@propagate_inbounds function add_bead_force!(force_energy::ForceEnergyFixedPoint{F_FBITS, E_FBITS}, i::Integer, force) where {F_FBITS, E_FBITS} + for j in 1:3 + force_energy.forces[3(i-1) + j] += convert_to_fixed_point(force[j], F_FBITS) + end + nothing +end +Base.@propagate_inbounds function zero_bead_force!(force_energy::ForceEnergyFixedPoint, i::Integer) + for j in 1:3 + force_energy.forces[3(i-1) + j] = 0 + end + nothing +end +function get_force!(force_energy::ForceEnergyFixedPoint{F_FBITS, E_FBITS}, force_out) where {F_FBITS, E_FBITS} + Base.require_one_based_indexing(force_out) + @argcheck 3*length(force_out) == length(force_energy.forces) + @inbounds for i in eachindex(force_out) + force_out[i] = @SVector [ + convert_to_floating_point(eltype(eltype(force_out)), force_energy.forces[3(i-1) + j], F_FBITS) + for j in 1:3 + ] + end + nothing +end +Base.@propagate_inbounds function get_bead_force(force_energy::ForceEnergyFixedPoint{F_FBITS, E_FBITS}, i::Integer, ::Type{F}=Float64)::SVector{3, F} where {F_FBITS, E_FBITS, F} + base = 3*(Int(i) - 1) + @SVector [ + convert_to_floating_point(F, force_energy.forces[base + j], F_FBITS) + for j in 1:3 + ] +end +function get_energy(force_energy::ForceEnergyFixedPoint{F_FBITS, E_FBITS}, ::Type{F}=Float64) where {F_FBITS, E_FBITS, F} + convert_to_floating_point(F, force_energy.energy, E_FBITS) +end +function combine_force_energy!( + fe_out::ForceEnergyFixedPoint{F, E}, + fe_in::ForceEnergyFixedPoint{F, E}, + ) where {F, E} + @argcheck length(fe_out.forces) == length(fe_in.forces) + fe_out.energy += fe_in.energy + @inbounds for i in eachindex(fe_in.forces) + fe_out.forces[i] += fe_in.forces[i] + end + fe_out +end + + +# Identical fixed-point force storage to `ForceEnergyFixedPoint`, but with no +# energy accumulator: `add_energy!` is a no-op and `get_energy` is intentionally +# not implemented, so reading the energy errors instead of returning a made-up +# value. +struct ForceNoEnergyFixedPoint{F_FBITS} <: ForceEnergy + forces::Memory{Int64} +end +function ForceNoEnergyFixedPoint{F_FBITS}(nbeads::Integer) where {F_FBITS} + forces = Memory{Int64}(undef, 3*nbeads) + fill!(forces, Int64(0)) + ForceNoEnergyFixedPoint{F_FBITS}(forces) +end +function zero_force_energy!(force_energy::ForceNoEnergyFixedPoint) + fill!(force_energy.forces, Int64(0)) + nothing +end +function get_nbeads(force_energy::ForceNoEnergyFixedPoint) + length(force_energy.forces)÷3 +end +function add_energy!(force_energy::ForceNoEnergyFixedPoint, energy) + nothing +end +function zero_energy!(force_energy::ForceNoEnergyFixedPoint) + nothing +end +Base.@propagate_inbounds function add_bead_force!(force_energy::ForceNoEnergyFixedPoint{F_FBITS}, i::Integer, force) where {F_FBITS} + for j in 1:3 + force_energy.forces[3(i-1) + j] += convert_to_fixed_point(force[j], F_FBITS) + end + nothing +end +Base.@propagate_inbounds function zero_bead_force!(force_energy::ForceNoEnergyFixedPoint, i::Integer) + for j in 1:3 + force_energy.forces[3(i-1) + j] = 0 + end + nothing +end +function get_force!(force_energy::ForceNoEnergyFixedPoint{F_FBITS}, force_out) where {F_FBITS} + Base.require_one_based_indexing(force_out) + @argcheck 3*length(force_out) == length(force_energy.forces) + @inbounds for i in eachindex(force_out) + force_out[i] = @SVector [ + convert_to_floating_point(eltype(eltype(force_out)), force_energy.forces[3(i-1) + j], F_FBITS) + for j in 1:3 + ] + end + nothing +end +Base.@propagate_inbounds function get_bead_force(force_energy::ForceNoEnergyFixedPoint{F_FBITS}, i::Integer, ::Type{F}=Float64)::SVector{3, F} where {F_FBITS, F} + base = 3*(Int(i) - 1) + @SVector [ + convert_to_floating_point(F, force_energy.forces[base + j], F_FBITS) + for j in 1:3 + ] +end +function combine_force_energy!( + fe_out::ForceNoEnergyFixedPoint{F}, + fe_in::ForceNoEnergyFixedPoint{F}, + ) where {F} + @argcheck length(fe_out.forces) == length(fe_in.forces) + @inbounds for i in eachindex(fe_in.forces) + fe_out.forces[i] += fe_in.forces[i] + end + fe_out +end + + +struct DebugForceEnergy{T} <: ForceEnergy + nbeads::Int64 + added_energy::Vector{T} + added_bead_force::Vector{Tuple{Int64, SVector{3, T}}} +end +function DebugForceEnergy{T}(nbeads::Integer) where T + DebugForceEnergy{T}(nbeads, [], []) +end +function zero_force_energy!(force_energy::DebugForceEnergy) + empty!(force_energy.added_energy) + empty!(force_energy.added_bead_force) + nothing +end +function get_nbeads(force_energy::DebugForceEnergy) + force_energy.nbeads +end +function add_energy!(force_energy::DebugForceEnergy, energy) + push!(force_energy.added_energy, energy) + nothing +end +function zero_energy!(force_energy::DebugForceEnergy) + empty!(force_energy.added_energy) + nothing +end +function add_bead_force!(force_energy::DebugForceEnergy, i::Integer, force) + checkbounds(1:force_energy.nbeads, i) + push!(force_energy.added_bead_force, (i, force)) + nothing +end +function zero_bead_force!(force_energy::DebugForceEnergy, i::Integer) + checkbounds(1:force_energy.nbeads, i) + filter!(((k, f),) -> k != i, force_energy.added_bead_force) + nothing +end +function get_force!(force_energy::DebugForceEnergy, force_out) + Base.require_one_based_indexing(force_out) + @argcheck length(force_out) == force_energy.nbeads + fill!(force_out, zero(eltype(force_out))) + for (i, f) in force_energy.added_bead_force + force_out[i] += f + end + nothing +end +function get_bead_force(force_energy::DebugForceEnergy{T}, i::Integer, ::Type{F}=Float64)::SVector{3, F} where {T, F} + checkbounds(1:force_energy.nbeads, i) + acc = zero(SVector{3, F}) + for (k, f) in force_energy.added_bead_force + if k == i + acc += SVector{3, F}(f) + end + end + acc +end +function get_energy(force_energy::DebugForceEnergy, ::Type{F}=Float64) where F + F(sum(force_energy.added_energy)) +end + + +# A fake force energy if no forces or energy need to be calculated. All writes +# are no-ops. Nothing is accumulated, so the getters (`get_force!`, +# `get_bead_force`, `get_energy`, ...) are intentionally not implemented and +# reading results errors instead of returning made-up values. +struct NullForceEnergy <: ForceEnergy + nbeads::Int64 +end +function zero_force_energy!(force_energy::NullForceEnergy) + nothing +end +function get_nbeads(force_energy::NullForceEnergy) + force_energy.nbeads +end +function add_energy!(force_energy::NullForceEnergy, energy) + nothing +end +function zero_energy!(force_energy::NullForceEnergy) + nothing +end +function add_bead_force!(force_energy::NullForceEnergy, i::Integer, force) + nothing +end +function zero_bead_force!(force_energy::NullForceEnergy, i::Integer) + nothing +end diff --git a/test/runtests.jl b/test/runtests.jl index 902c4a8..1c171a9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,3 +8,5 @@ include("test-pointcelllist.jl") include("test-linesegcelllist.jl") +include("test-force-energy.jl") + diff --git a/test/test-force-energy.jl b/test/test-force-energy.jl new file mode 100644 index 0000000..a99a6ac --- /dev/null +++ b/test/test-force-energy.jl @@ -0,0 +1,292 @@ +using SimplexCellLists +using Test +using StaticArrays + +@testset "ForceEnergy" begin + @testset "ForceEnergyT{Float64}" begin + nbeads = 4 + fe = ForceEnergyFloat64(nbeads) + + # Initially zero on every bead. + for i in 1:nbeads + @test get_bead_force(fe, i) === SA[0.0, 0.0, 0.0] + end + + # Single push picks up on get_bead_force. + add_bead_force!(fe, 2, SA[1.5, -2.0, 3.25]) + @test get_bead_force(fe, 1) === SA[0.0, 0.0, 0.0] + @test get_bead_force(fe, 2) === SA[1.5, -2.0, 3.25] + @test get_bead_force(fe, 3) === SA[0.0, 0.0, 0.0] + + # Accumulates across multiple pushes to the same bead. + add_bead_force!(fe, 2, SA[0.5, 1.0, -0.25]) + @test get_bead_force(fe, 2) === SA[2.0, -1.0, 3.0] + + # Matches get_force! output for the same index. + out = zeros(SVector{3, Float64}, nbeads) + get_force!(fe, out) + for i in 1:nbeads + @test get_bead_force(fe, i) === out[i] + end + + # Energy accumulates, and zero_energy! resets only the energy. + add_energy!(fe, 1.5) + add_energy!(fe, -0.25) + @test get_energy(fe) === 1.25 + zero_energy!(fe) + @test get_energy(fe) === 0.0 + @test get_bead_force(fe, 2) === SA[2.0, -1.0, 3.0] + + # zero_bead_force! wipes only that bead. + add_bead_force!(fe, 3, SA[1.0, 1.0, 1.0]) + zero_bead_force!(fe, 2) + @test get_bead_force(fe, 2) === SA[0.0, 0.0, 0.0] + @test get_bead_force(fe, 3) === SA[1.0, 1.0, 1.0] + + # zero_force_energy! wipes everything. + add_energy!(fe, 1.0) + zero_force_energy!(fe) + @test get_energy(fe) === 0.0 + for i in 1:nbeads + @test get_bead_force(fe, i) === SA[0.0, 0.0, 0.0] + end + end + + @testset "ForceEnergyT{Float32}" begin + # Returned type is Float64 even when the backing storage is Float32. + fe = ForceEnergyFloat32(2) + add_bead_force!(fe, 1, SA{Float32}[1.0f0, 2.0f0, -3.0f0]) + f = get_bead_force(fe, 1) + @test f isa SVector{3, Float64} + @test f == SA[1.0, 2.0, -3.0] + end + + @testset "ForceEnergyFixedPoint" begin + nbeads = 3 + fe = ForceEnergyFixedPoint{30, 30}(nbeads) + + for i in 1:nbeads + @test get_bead_force(fe, i) === SA[0.0, 0.0, 0.0] + end + + add_bead_force!(fe, 1, SA[0.25, -0.5, 1.0]) + add_bead_force!(fe, 3, SA[2.0, 0.0, -1.5]) + add_bead_force!(fe, 1, SA[0.125, 0.25, 0.0]) + + f1 = get_bead_force(fe, 1) + @test f1 isa SVector{3, Float64} + @test f1 ≈ SA[0.375, -0.25, 1.0] + @test get_bead_force(fe, 2) === SA[0.0, 0.0, 0.0] + @test get_bead_force(fe, 3) ≈ SA[2.0, 0.0, -1.5] + + # Matches get_force! output. + out = zeros(SVector{3, Float64}, nbeads) + get_force!(fe, out) + for i in 1:nbeads + @test get_bead_force(fe, i) ≈ out[i] + end + + # Energy accumulates in fixed point, and zero_energy! resets it. + add_energy!(fe, 1.5) + add_energy!(fe, -0.25) + @test get_energy(fe) ≈ 1.25 + @test get_energy(fe, Float32) isa Float32 + zero_energy!(fe) + @test get_energy(fe) === 0.0 + + # zero_bead_force! wipes only that bead. + zero_bead_force!(fe, 1) + @test get_bead_force(fe, 1) === SA[0.0, 0.0, 0.0] + @test get_bead_force(fe, 3) ≈ SA[2.0, 0.0, -1.5] + + # Zeroing wipes per-bead view too. + add_energy!(fe, 1.0) + zero_force_energy!(fe) + @test get_energy(fe) === 0.0 + for i in 1:nbeads + @test get_bead_force(fe, i) === SA[0.0, 0.0, 0.0] + end + end + + @testset "ForceNoEnergyFixedPoint" begin + nbeads = 3 + fe = ForceNoEnergyFixedPoint{30}(nbeads) + + @test get_nbeads(fe) == nbeads + for i in 1:nbeads + @test get_bead_force(fe, i) === SA[0.0, 0.0, 0.0] + end + + add_bead_force!(fe, 1, SA[0.25, -0.5, 1.0]) + add_bead_force!(fe, 3, SA[2.0, 0.0, -1.5]) + add_bead_force!(fe, 1, SA[0.125, 0.25, 0.0]) + @test get_bead_force(fe, 1) ≈ SA[0.375, -0.25, 1.0] + @test get_bead_force(fe, 2) === SA[0.0, 0.0, 0.0] + @test get_bead_force(fe, 3) ≈ SA[2.0, 0.0, -1.5] + + # Energy is discarded: add_energy!/zero_energy! are no-ops, and reading + # the energy errors instead of returning a made-up value. + add_energy!(fe, 5.0) + zero_energy!(fe) + @test_throws MethodError get_energy(fe) + + # Matches get_force! output. + out = zeros(SVector{3, Float64}, nbeads) + get_force!(fe, out) + for i in 1:nbeads + @test get_bead_force(fe, i) ≈ out[i] + end + + # zero_bead_force! wipes only that bead. + zero_bead_force!(fe, 1) + @test get_bead_force(fe, 1) === SA[0.0, 0.0, 0.0] + @test get_bead_force(fe, 3) ≈ SA[2.0, 0.0, -1.5] + + zero_force_energy!(fe) + for i in 1:nbeads + @test get_bead_force(fe, i) === SA[0.0, 0.0, 0.0] + end + end + + @testset "combine_force_energy! bit-exact for matching ForceNoEnergyFixedPoint" begin + nbeads = 4 + a = ForceNoEnergyFixedPoint{30}(nbeads) + b = ForceNoEnergyFixedPoint{30}(nbeads) + ref = ForceNoEnergyFixedPoint{30}(nbeads) + + contribs_a = [(1, SA[0.25, -0.5, 1.0]), (3, SA[2.0, 0.0, -1.5]), (1, SA[0.125, 0.25, 0.0])] + contribs_b = [(2, SA[0.5, 0.5, 0.5]), (3, SA[-1.0, 1.0, 0.0]), (4, SA[0.75, -0.25, 0.125])] + + for (i, f) in contribs_a + add_bead_force!(a, i, f) + add_bead_force!(ref, i, f) + end + for (i, f) in contribs_b + add_bead_force!(b, i, f) + add_bead_force!(ref, i, f) + end + + combine_force_energy!(a, b) + + # Underlying Int64 storage matches bit-for-bit. + @test a.forces == ref.forces + end + + @testset "combine_force_energy! bit-exact for matching FixedPoint" begin + # Combining two ForceEnergyFixedPoint buffers with matching F,E should + # produce results identical (bit-for-bit) to accumulating every + # contribution into a single buffer from the start. + nbeads = 4 + F, E = 30, 30 + a = ForceEnergyFixedPoint{F, E}(nbeads) + b = ForceEnergyFixedPoint{F, E}(nbeads) + ref = ForceEnergyFixedPoint{F, E}(nbeads) + + contribs_a = [(1, SA[0.25, -0.5, 1.0]), (3, SA[2.0, 0.0, -1.5]), (1, SA[0.125, 0.25, 0.0])] + contribs_b = [(2, SA[0.5, 0.5, 0.5]), (3, SA[-1.0, 1.0, 0.0]), (4, SA[0.75, -0.25, 0.125])] + + for (i, f) in contribs_a + add_bead_force!(a, i, f) + add_bead_force!(ref, i, f) + end + add_energy!(a, 1.5) + add_energy!(ref, 1.5) + + for (i, f) in contribs_b + add_bead_force!(b, i, f) + add_bead_force!(ref, i, f) + end + add_energy!(b, -0.75) + add_energy!(ref, -0.75) + + combine_force_energy!(a, b) + + # Underlying Int64 storage matches bit-for-bit. + @test a.forces == ref.forces + @test a.energy[] === ref.energy[] + end + + @testset "combine_force_energy! generic path on FixedPoint" begin + # Using the generic combine_force_energy! (e.g. via differing F,E or a + # mixed-type combine) goes through float round-tripping. Sanity-check + # that the result still agrees on the values we use here, but is *not* + # guaranteed bit-exact. + nbeads = 2 + a = ForceEnergyFixedPoint{30, 30}(nbeads) + b = ForceEnergyFixedPoint{20, 30}(nbeads) + add_bead_force!(a, 1, SA[0.25, 0.5, -1.0]) + add_bead_force!(b, 1, SA[0.125, -0.25, 2.0]) + combine_force_energy!(a, b) + @test get_bead_force(a, 1) ≈ SA[0.375, 0.25, 1.0] + end + + @testset "DebugForceEnergy" begin + nbeads = 3 + fe = DebugForceEnergy{Float64}(nbeads) + + for i in 1:nbeads + @test get_bead_force(fe, i) === SA[0.0, 0.0, 0.0] + end + + add_bead_force!(fe, 2, SA[1.0, 2.0, 3.0]) + add_bead_force!(fe, 2, SA[0.5, -1.0, 0.0]) + add_bead_force!(fe, 3, SA[-1.0, 0.0, 0.0]) + + @test get_bead_force(fe, 1) === SA[0.0, 0.0, 0.0] + @test get_bead_force(fe, 2) === SA[1.5, 1.0, 3.0] + @test get_bead_force(fe, 3) === SA[-1.0, 0.0, 0.0] + + # Matches get_force! output. + out = zeros(SVector{3, Float64}, nbeads) + get_force!(fe, out) + for i in 1:nbeads + @test get_bead_force(fe, i) === out[i] + end + + # Energy sums the recorded contributions, zero_energy! clears them. + add_energy!(fe, 1.5) + add_energy!(fe, -0.25) + @test get_energy(fe) === 1.25 + zero_energy!(fe) + @test get_energy(fe) === 0.0 + + # zero_bead_force! wipes only that bead. + zero_bead_force!(fe, 2) + @test get_bead_force(fe, 2) === SA[0.0, 0.0, 0.0] + @test get_bead_force(fe, 3) === SA[-1.0, 0.0, 0.0] + + # Out-of-range index is rejected. + @test_throws Exception get_bead_force(fe, 0) + @test_throws Exception get_bead_force(fe, nbeads + 1) + @test_throws Exception add_bead_force!(fe, nbeads + 1, SA[1.0, 0.0, 0.0]) + @test_throws Exception zero_bead_force!(fe, nbeads + 1) + + add_energy!(fe, 1.0) + zero_force_energy!(fe) + @test get_energy(fe) === 0.0 + for i in 1:nbeads + @test get_bead_force(fe, i) === SA[0.0, 0.0, 0.0] + end + end + + @testset "NullForceEnergy" begin + nbeads = 4 + fe = NullForceEnergy(nbeads) + + @test get_nbeads(fe) == nbeads + + # All writes are no-ops. + add_bead_force!(fe, 2, SA[1.0, 2.0, 3.0]) + add_energy!(fe, 1.5) + zero_bead_force!(fe, 2) + zero_energy!(fe) + zero_force_energy!(fe) + + # Nothing is accumulated, so reading results errors. + @test_throws MethodError get_energy(fe) + @test_throws MethodError get_bead_force(fe, 1) + @test_throws MethodError get_force!(fe, zeros(SVector{3, Float64}, nbeads)) + @test_throws MethodError get_force(fe) + end +end +nothing diff --git a/test/test-mindistance.jl b/test/test-mindistance.jl index c927374..76eab9d 100644 --- a/test/test-mindistance.jl +++ b/test/test-mindistance.jl @@ -75,3 +75,4 @@ end # This issue was resolved in https://github.com/oxfordcontrol/Clarabel.jl/pull/114 @test ref_d2min ≈ d2min atol = 1E-11 rtol = 1E-11 end +nothing From de5deb4e3dbfab0d098645ac6fc30eac57474f03 Mon Sep 17 00:00:00 2001 From: nhz2 Date: Mon, 10 Aug 2026 10:20:11 +0200 Subject: [PATCH 2/3] add docstrings --- src/SimplexCellLists.jl | 5 +- src/force-energy.jl | 149 +++++++++++++++++++++++++++++++++++--- test/test-force-energy.jl | 13 +++- 3 files changed, 154 insertions(+), 13 deletions(-) diff --git a/src/SimplexCellLists.jl b/src/SimplexCellLists.jl index ce85b80..8648510 100644 --- a/src/SimplexCellLists.jl +++ b/src/SimplexCellLists.jl @@ -47,6 +47,7 @@ let T = Int64, F = Float32 end include("force-energy.jl") +export ForceEnergy export ForceEnergyT export ForceEnergyFloat64 export ForceEnergyFloat32 @@ -66,6 +67,4 @@ export get_energy export get_force export combine_force_energy! - - -end \ No newline at end of file +end diff --git a/src/force-energy.jl b/src/force-energy.jl index 248906f..f7baa3f 100644 --- a/src/force-energy.jl +++ b/src/force-energy.jl @@ -38,11 +38,89 @@ size, e.g. to reduce per-thread partials). """ abstract type ForceEnergy end -function get_force(force_energy::ForceEnergy, F=Float64) +""" + add_bead_force!(force_energy::ForceEnergy, i::Integer, force)::Nothing + +Add `force`, a 3-element vector, to the accumulated force of bead `i`. +""" +function add_bead_force! end + +""" + zero_bead_force!(force_energy::ForceEnergy, i::Integer)::Nothing + +Set the accumulated force of bead `i` to zero. +""" +function zero_bead_force! end + +""" + add_energy!(force_energy::ForceEnergy, energy)::Nothing + +Add the scalar `energy` to the accumulated total energy. +""" +function add_energy! end + +""" + zero_energy!(force_energy::ForceEnergy)::Nothing + +Set the accumulated total energy to zero, leaving forces unchanged. +""" +function zero_energy! end + +""" + get_nbeads(force_energy::ForceEnergy)::Int + +Return the number of beads. +""" +function get_nbeads end + +""" + zero_force_energy!(force_energy::ForceEnergy)::Nothing + +Set all accumulated forces and the total energy to zero. +""" +function zero_force_energy! end + +""" + get_force!(force_energy::ForceEnergy, force_out)::Nothing + +Write the accumulated per-bead forces into `force_out`, a one-based vector of +3-element vectors of length [`get_nbeads`](@ref)`(force_energy)`. +""" +function get_force! end + +""" + get_bead_force(force_energy::ForceEnergy, i::Integer, F::Type=Float64)::SVector{3, F} + +Return the accumulated force of bead `i`, converted to element type `F`. +""" +function get_bead_force end + +""" + get_energy(force_energy::ForceEnergy, F::Type=Float64)::F + +Return the accumulated total energy, converted to type `F`. +""" +function get_energy end + +""" + get_force(force_energy::ForceEnergy, F::Type=Float64)::Vector{SVector{3, F}} + +Return a new vector of the accumulated per-bead forces. +""" +function get_force(force_energy::ForceEnergy, ::Type{F}=Float64) where F force_out = zeros(SVector{3, F}, get_nbeads(force_energy)) get_force!(force_energy, force_out) force_out end + +""" + combine_force_energy!(fe_out::ForceEnergy, fe_in::ForceEnergy) -> fe_out + +Merge `fe_in` into `fe_out`, which must have the same number of beads, by +adding its accumulated forces and energy. Useful to reduce per-thread partial +accumulators. The generic fallback round-trips through `Float64`; matching +fixed-point types merge bit-exactly. +""" function combine_force_energy!(fe_out::ForceEnergy, fe_in::ForceEnergy) @argcheck get_nbeads(fe_in) == get_nbeads(fe_out) add_energy!(fe_out, get_energy(fe_in)) @@ -52,11 +130,29 @@ function combine_force_energy!(fe_out::ForceEnergy, fe_in::ForceEnergy) fe_out end +""" + ForceEnergyT{T}(nbeads::Integer) + +A [`ForceEnergy`](@ref) accumulator storing per-bead forces and the total +energy as numbers of type `T`, summed in place. If `T` is an inexact type like +`Float64`, results can depend on the order contributions are added in, so +parallel reductions are not reproducible bit-for-bit. +""" mutable struct ForceEnergyT{T} <: ForceEnergy const forces::Memory{SVector{3,T}} energy::T end +""" + ForceEnergyFloat64 + +Alias for [`ForceEnergyT`](@ref)`{Float64}`. +""" const ForceEnergyFloat64 = ForceEnergyT{Float64} +""" + ForceEnergyFloat32 + +Alias for [`ForceEnergyT`](@ref)`{Float32}`. +""" const ForceEnergyFloat32 = ForceEnergyT{Float32} function ForceEnergyT{T}(nbeads::Integer) where T forces = Memory{SVector{3,T}}(undef, nbeads) @@ -106,6 +202,23 @@ end @inline convert_to_fixed_point(x, fbits) = unsafe_trunc(Int64, x*(oftype(x, 2)^fbits)) @inline convert_to_floating_point(T, x, fbits) = T(x)*(convert(T, 2)^-fbits) +""" + ForceEnergyFixedPoint{F_FBITS, E_FBITS}(nbeads::Integer) + +A [`ForceEnergy`](@ref) accumulator storing per-bead forces and the total +energy as `Int64` fixed-point numbers with `F_FBITS` and `E_FBITS` fractional +bits respectively. Accumulation is exact once inputs are quantized, so results +are independent of the order contributions are added in, and accumulators with +matching parameters merge bit-exactly with [`combine_force_energy!`](@ref) +(e.g. to reduce per-thread partials reproducibly). + +Added values are quantized by truncation toward zero, losing up to +`2^-F_FBITS` (or `2^-E_FBITS`) of magnitude per contribution. Negating an +input exactly negates its quantized value, so equal-and-opposite contributions +cancel exactly. Inputs must be finite with magnitude less than +`2^(63 - FBITS)`; non-finite or out-of-range values silently corrupt the +accumulator (`unsafe_trunc`), as does overflow of the running `Int64` sums. +""" mutable struct ForceEnergyFixedPoint{F_FBITS, E_FBITS} <: ForceEnergy const forces::Memory{Int64} energy::Int64 @@ -177,10 +290,14 @@ function combine_force_energy!( end -# Identical fixed-point force storage to `ForceEnergyFixedPoint`, but with no -# energy accumulator: `add_energy!` is a no-op and `get_energy` is intentionally -# not implemented, so reading the energy errors instead of returning a made-up -# value. +""" + ForceNoEnergyFixedPoint{F_FBITS}(nbeads::Integer) + +Identical fixed-point force storage to [`ForceEnergyFixedPoint`](@ref), but +with no energy accumulator: [`add_energy!`](@ref) is a no-op and +[`get_energy`](@ref) is intentionally not implemented, so reading the energy +errors instead of returning a made-up value. +""" struct ForceNoEnergyFixedPoint{F_FBITS} <: ForceEnergy forces::Memory{Int64} end @@ -244,6 +361,15 @@ function combine_force_energy!( end +""" + DebugForceEnergy{T}(nbeads::Integer) + +A [`ForceEnergy`](@ref) accumulator for testing that records every +[`add_energy!`](@ref) and [`add_bead_force!`](@ref) contribution individually +as type `T` instead of summing in place. Readers sum the recorded +contributions on demand, and [`zero_bead_force!`](@ref) removes the recorded +contributions to that bead. All bead indices are bounds-checked. +""" struct DebugForceEnergy{T} <: ForceEnergy nbeads::Int64 added_energy::Vector{T} @@ -302,10 +428,15 @@ function get_energy(force_energy::DebugForceEnergy, ::Type{F}=Float64) where F end -# A fake force energy if no forces or energy need to be calculated. All writes -# are no-ops. Nothing is accumulated, so the getters (`get_force!`, -# `get_bead_force`, `get_energy`, ...) are intentionally not implemented and -# reading results errors instead of returning made-up values. +""" + NullForceEnergy(nbeads::Integer) + +A fake [`ForceEnergy`](@ref) if no forces or energy need to be calculated. All +writes are no-ops. Nothing is accumulated, so the getters +([`get_force!`](@ref), [`get_bead_force`](@ref), [`get_energy`](@ref), ...) +are intentionally not implemented and reading results errors instead of +returning made-up values. +""" struct NullForceEnergy <: ForceEnergy nbeads::Int64 end diff --git a/test/test-force-energy.jl b/test/test-force-energy.jl index a99a6ac..afd49aa 100644 --- a/test/test-force-energy.jl +++ b/test/test-force-energy.jl @@ -106,6 +106,17 @@ using StaticArrays for i in 1:nbeads @test get_bead_force(fe, i) === SA[0.0, 0.0, 0.0] end + + # Quantization truncates toward zero symmetrically: adding a + # contribution and its exact negation cancels exactly, even for values + # that don't quantize exactly. + v = SA[0.3, -1e-9, 7.123456789] + add_bead_force!(fe, 1, v) + add_bead_force!(fe, 1, -v) + add_energy!(fe, 0.3) + add_energy!(fe, -0.3) + @test get_bead_force(fe, 1) === SA[0.0, 0.0, 0.0] + @test get_energy(fe) === 0.0 end @testset "ForceNoEnergyFixedPoint" begin @@ -203,7 +214,7 @@ using StaticArrays # Underlying Int64 storage matches bit-for-bit. @test a.forces == ref.forces - @test a.energy[] === ref.energy[] + @test a.energy === ref.energy end @testset "combine_force_energy! generic path on FixedPoint" begin From c236b4ea07ff60ce899e7c38dbe32f1d82503939 Mon Sep 17 00:00:00 2001 From: nhz2 Date: Mon, 10 Aug 2026 12:10:24 +0200 Subject: [PATCH 3/3] precompile and more test coverage --- src/SimplexCellLists.jl | 22 ++++++++++++++++++++++ src/force-energy.jl | 2 +- test/test-force-energy.jl | 36 ++++++++++++++++++++++++++++++++++++ 3 files changed, 59 insertions(+), 1 deletion(-) diff --git a/src/SimplexCellLists.jl b/src/SimplexCellLists.jl index 8648510..adf253e 100644 --- a/src/SimplexCellLists.jl +++ b/src/SimplexCellLists.jl @@ -67,4 +67,26 @@ export get_energy export get_force export combine_force_energy! +# Other fixed point fractional bit choices must be specialized on the +# caller's parameters. +for (FE, T) in ( + (ForceEnergyFloat64, Float64), + (ForceEnergyFixedPoint{30, 30}, Float64), + (ForceNoEnergyFixedPoint{30}, Float64), + ) + precompile(FE, (Int,)) + precompile(zero_force_energy!, (FE,)) + precompile(add_bead_force!, (FE, Int, SVector{3,T})) + precompile(zero_bead_force!, (FE, Int)) + precompile(add_energy!, (FE, T)) + precompile(zero_energy!, (FE,)) + precompile(get_force!, (FE, Vector{SVector{3,T}})) + precompile(get_bead_force, (FE, Int)) + precompile(get_force, (FE,)) + precompile(combine_force_energy!, (FE, FE)) +end +# ForceNoEnergyFixedPoint has no get_energy method. +precompile(get_energy, (ForceEnergyFloat64,)) +precompile(get_energy, (ForceEnergyFixedPoint{30, 30},)) + end diff --git a/src/force-energy.jl b/src/force-energy.jl index f7baa3f..2b7c710 100644 --- a/src/force-energy.jl +++ b/src/force-energy.jl @@ -376,7 +376,7 @@ struct DebugForceEnergy{T} <: ForceEnergy added_bead_force::Vector{Tuple{Int64, SVector{3, T}}} end function DebugForceEnergy{T}(nbeads::Integer) where T - DebugForceEnergy{T}(nbeads, [], []) + DebugForceEnergy{T}(nbeads, T[], Tuple{Int64, SVector{3, T}}[]) end function zero_force_energy!(force_energy::DebugForceEnergy) empty!(force_energy.added_energy) diff --git a/test/test-force-energy.jl b/test/test-force-energy.jl index afd49aa..d761792 100644 --- a/test/test-force-energy.jl +++ b/test/test-force-energy.jl @@ -29,6 +29,11 @@ using StaticArrays @test get_bead_force(fe, i) === out[i] end + # Allocating get_force matches, in the requested element type. + @test get_force(fe) == out + @test get_force(fe) isa Vector{SVector{3, Float64}} + @test get_force(fe, Float32) isa Vector{SVector{3, Float32}} + # Energy accumulates, and zero_energy! resets only the energy. add_energy!(fe, 1.5) add_energy!(fe, -0.25) @@ -86,6 +91,13 @@ using StaticArrays @test get_bead_force(fe, i) ≈ out[i] end + # get_force! converts to the output vector's element type. + out32 = zeros(SVector{3, Float32}, nbeads) + get_force!(fe, out32) + for i in 1:nbeads + @test out32[i] ≈ SVector{3, Float32}(get_bead_force(fe, i)) + end + # Energy accumulates in fixed point, and zero_energy! resets it. add_energy!(fe, 1.5) add_energy!(fe, -0.25) @@ -231,6 +243,30 @@ using StaticArrays @test get_bead_force(a, 1) ≈ SA[0.375, 0.25, 1.0] end + @testset "combine_force_energy! generic path on ForceEnergyT" begin + nbeads = 3 + a = ForceEnergyFloat64(nbeads) + b = ForceEnergyFloat64(nbeads) + add_bead_force!(a, 1, SA[1.0, 2.0, 3.0]) + add_energy!(a, 0.5) + add_bead_force!(b, 1, SA[0.5, -1.0, 0.25]) + add_bead_force!(b, 2, SA[1.0, 0.0, 0.0]) + add_energy!(b, 0.25) + + @test combine_force_energy!(a, b) === a + @test get_bead_force(a, 1) === SA[1.5, 1.0, 3.25] + @test get_bead_force(a, 2) === SA[1.0, 0.0, 0.0] + @test get_bead_force(a, 3) === SA[0.0, 0.0, 0.0] + @test get_energy(a) === 0.75 + + # The input accumulator is unchanged. + @test get_bead_force(b, 1) === SA[0.5, -1.0, 0.25] + @test get_energy(b) === 0.25 + + # Mismatched bead counts are rejected. + @test_throws Exception combine_force_energy!(a, ForceEnergyFloat64(nbeads + 1)) + end + @testset "DebugForceEnergy" begin nbeads = 3 fe = DebugForceEnergy{Float64}(nbeads)