Repository navigation
BREAKING: SIMD collide forces - #27
Merged
Merged
Conversation
Policies now define force laws with `nl_edge_force`, which takes the two simplexes, relative to the first vertex of `x` and converted to `calc_type`, plus their `dist_sqr_r_s_t` result, and returns the energy and the forces on each vertex. It is pure math with no memory access, so it can be vectorized or called from a GPU kernel. `nl_edge_forces!` loads the positions, does the `d2 < L^2` range check, calls the hook and scatters the forces. - Point-point, point-line and point-triangle share one force law, spread over the vertices by the barycentric weights. - Line-line keeps the switchover blend with the boundary critical points. `line_line_force_core` is removed. - `calc_type` is now a type, like `Float32`, and other types such as `BigFloat` use the same scalar path. - `nl_edge_force` is public. - The version goes to 0.3.0-dev, since this breaks the API. Add tests of Float32 forces for large simplexes 0.01 apart, and of the scalar path with `BigFloat`. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
For `calc_type` `Float32` or `Float64`, `nl_forces!` checks the range
8 edges at a time, with `dist_sqr` on `SIMD.Vec` lanes. Each edge in
range is copied to an `EdgeBuffer`, a structure of arrays in `Memory`,
and every 8 edges the forces are computed in a plain loop over the
lanes, which LLVM vectorizes, then added to the accumulator in edge
order. The last edges, fewer than 8, take the scalar path. The results
are the same as the scalar path, as both call the same functions on
the same values.
- The positions, `L`, and the forces share one `Memory{T}`, with the
forces stored over the positions. With separate buffers LLVM needs
runtime alias checks and leaves the force loop scalar.
- The params and index parts get their own `Memory`, so params need
not be isbits. The params are loaded to local variables before the
force loop.
- The buffers stay under about 2000 bytes, so Julia keeps them on the
stack and `collide_forces!` doesn't allocate.
The benchmark now uses positions in the calculation type. Against
main, collide_forces! is 4.4 to 5.3 times faster with Float64 and 3.9
to 4.7 times faster with Float32.
Add tests that the vectorized path matches the scalar path exactly,
doesn't allocate, and works with params that aren't isbits.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
A zero length line segment at a distance gave NaN forces with a finite energy: the switchover derivative divided by `p*q`, `p` and `q` in turn, which are zero. Now there is one inverse, `ipq = 1/(p*q)`, which is zero when `p*q` is at most `floatmin`. Then the switchover is zero, the limit as the line segment shrinks, so the energy is continuous, and the switchover gives no force. `q/p + b*g/a` is written as `c*ipq*(1 + g^2)`, and `p/q + b*g/c` as `a*ipq*(1 + g^2)`. `Δ` is now `‖P×Q‖^2` from `cross3`, the `n` of `dist_sqr_r_s_t`, instead of `a*c - b*b`. That cancelled for parallel line segments and could go negative, which made the switchover slightly negative. The forces agree with the old ones to about 4e-15 elsewhere, and match finite differences down to parallel line segments. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
When the same positions are run over and over, the branch predictor learns which groups of 8 edges have one in range, if there are fewer than about 6000 such groups. The ring has about 300, so its time was too low: hits had almost no mispredicts, against about 1.3 per group of edges with a hit otherwise. It now runs 32 copies jittered by 0.3 nm in turn, one call per sample, which takes it from 20.4 to 22.4 µs with Float64. The isosphere systems have enough groups already, so they keep one copy. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Compute the edge parameters with divisions and the separations with fma3, as dist_sqr_r_s_t does, so the end point candidates of x match its r0 and r1 exactly. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The comment said the switchover gives no force in the limit as a line segment shrinks. The energy has that limit, but the gradient of the switchover along the short line segment stays about `‖Q‖*sin²θ*switchover_scale`, so the force depends on its direction. Zero is a choice, not the limit. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
`cross3(P, P)` is a rounding residual of about `eps*p*q`, not zero, which switched over parallel diagonal line segments in Float32 and gave about 21% too little force. Use `cross3_accurate`, which cancels it exactly. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
It fails under CI's --check-bounds=yes setting, and the benchmarks already check performance. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The expanded gradient `c/(p*q)*(1 + g^2)*P - 2*g*Q` cancels to a rounding of `‖Q‖` for parallel line segments, so it gave a force along them even when the accurate `P×Q` was exactly zero: 3.8e-6 in Float32 for a pair pushed apart by 1e-6. Writing it as `(2*Q×n - Δ/a*P)/(p*q)`, and the same for `Q`, makes both terms scale with `‖n‖`, so the gradient is exactly zero for parallel line segments. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This PR speeds up collide forces using SIMD. This is a breaking change for custom force laws. The inner force function needs to be changed from
nl_edge_forces!tonl_edge_forces.