Objective

A big step in #319.
Replaces #706.

Currently, Avian's constraint solver is entirely single-threaded and uses scalar operations. A naive par_iter approach is not possible with a Gauss-Seidel solver without race conditions.

One form of solver parallelism is inter-island parallelism, meaning that bodies are assigned to islands based on the constraint graph, and each island can be safely solved in parallel, since they are guaranteed to be fully disjoint. Rapier uses this for its solver.

However, inter-island parallelism does not help with large islands, such as large piles of bodies. Constraints within each island must still be solved serially. Additionally, it does not help in vectorizing the constraint solver via wide SIMD. Simulation islands are therefore not enough for the performance we are aiming for (though we do still need them for sleeping and waking!).

State-of-the-art physics solvers instead use various forms of graph coloring. The idea is that each constraint is assigned a color such that no adjacent edge shares the same color. This ensures that each body only appears within each color once at most, and all constraints within a given color can be safely solved in parallel without race conditions. This not only makes it trivial to solve constraints with multi-threading, but also opens the door to solving constraints with wide SIMD.

The scope of this PR is to implement a graph coloring algorithm for contact constraints, and to make any related important changes around contact management. Wide SIMD and joint graph coloring are left to follow-ups.

Solution

• Rework the ContactGraph to work better with graph coloring and optimize contact management.
  • The graph is now a StableUnGraph that stores ContactEdges with stable ContactIds. Contact edges store cold connectivity data.
  • Store the actual contact pairs in separate active_pairs and sleeping_pairs vectors. Only active pairs are considered in graph coloring and iterated in the narrow phase.
• Move ContactConstraint preparation out of the narrow phase and into the solver.
• Implement a ConstraintGraph for greedy graph coloring.
  • Each ContactManifold has its own constraint, so the coloring is done at the manifold-level.
  • Each GraphColor stores ContactManifoldHandles to map to contact manifolds in the ContactGraph.
  • Each ContactEdge stores ContactConstraintHandles to map to contact constraints in the ConstraintGraph.
  • Graph colors use BitVecs to track which bodies are assigned to them.

References

The original design was largely based on Box2D's graph coloring implementation. However, unlike Box2D, we can have multiple contact manifolds per contact pair (for e.g. triangle meshes or other composite shapes), so we need to perform coloring at the manifold-level, which complicates things a bit. That part is my own design. Additionally, we store contact pairs in the contact graph instead of the constraint graph.

Some useful references:

• Graph coloring
• High-Performance Physical Simulations on Next-Generation Architecture with Many Cores
• Box2D - SIMD Matters by Erin Catto

Performance

Running the Large Pyramid 2D benchmark with this branch, the main branch, and bevy_rapier2d 0.30 (configured to match Avian), we get the following results:

This means that with graph coloring, Avian is now faster than Rapier for this scene! Rapier's profile is very flat, as it does not have solver parallelism for large islands.

Note that Rapier is still generally faster than Avian for scenes that actually benefit from inter-island parallelism. For example, the Many Pyramids 2D benchmark looks like this:

However, judging by the shapes of the graphs, our multi-threaded scaling still seems better. I think Avian is slower here not because of the solver, but because of other less optimized parts, like the broad phase and spatial query pipeline updates. The upcoming BVH broad phase and spatial query rework should largely resolve this.

Here you can see a more detailed breakdown of Avian's performance on the main branch:

And now with graph coloring, on this branch:

At the top of the profile you can see the contacts being assigned to graph colors. You can see how the narrow phase is now significantly faster, and more importantly the solver is over 3x as fast on my machine, resulting in nearly 2x total performance compared to main. This difference tends to grow the more contacts there are.

Something that is not shown in this profile is how contacts are now split into separate sets for active and sleeping bodies. It means that the narrow phase can entirely skip iteration over sleeping contacts, which may be a noticeable improvement for larger scenes with lots of sleeping bodies.

Future Work

• Wide SIMD for contacts
• Graph coloring for joints
• Simulation islands for sleeping and waking
• Split sensor logic to be separate from the narrow phase

Migration Guide

Contact graph

• The ContactGraph now stores ContactEdges in the graph, containing only edge connectivity information with stable ContactIds.
• The ContactGraph now stores ContactPairs in separate active_pairs and sleeping_pairs lists.
• iter and iter_mut have been renamed to iter_active and iter_active_mut / iter_sleeping and iter_sleeping_mut.
• iter_touching and iter_touching_mut have been renamed to iter_touching_active and iter_touching_active_mut / iter_touching_sleeping and iter_touching_sleeping_mut.
• collisions_with has been renamed to contact_pairs_with.
• collisions_with_mut has been removed for now.
• add_pair has been replaced by add_edge, and add_pair_with_key has been replaced by add_edge_and_key_with.
• insert_pair and insert_pair_with have been removed. Add or update pairs manually instead.
• remove_pair has been renamed to remove_edge.

Contact types

• ContactPair now stores the ContactId associated with the ContactEdge in the ContactGraph. It must be provided to ContactPair::new. It is initialized automatically by calls to methods like ContactGraph::add_edge.
• ContactPairFlags no longer stores whether CONTACT_EVENTS are enabled. It is instead stored on ContactEdgeFlags.
• ContactManifold no longer stores the index of the manifold.

Determinism

The ContactGraph uses a StableUnGraph with stable indices. A side effect of this is that the order of old contacts can affect where new contacts are added. Thus, if you reset a scene without clearing the contact graph, you may notice slightly different behavior.

For consistent, deterministic behavior when resetting scenes, consider calling Collisions::clear to also reset the contact graph and make sure physics starts with a clean slate.