# 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](https://box2d.org/posts/2023/10/simulation-islands/)*, 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](https://en.wikipedia.org/wiki/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 `ContactEdge`s with stable `ContactId`s. 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 `ContactManifoldHandle`s to map to contact manifolds in the `ContactGraph`.
  - Each `ContactEdge` stores `ContactConstraintHandle`s to map to contact constraints in the `ConstraintGraph`.
  - Graph colors use `BitVec`s 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](https://en.wikipedia.org/wiki/Graph_coloring)
- [High-Performance Physical Simulations on Next-Generation Architecture with Many Cores][Intel paper]
- [Box2D - SIMD Matters] by [Erin Catto]

[Intel paper]: http://web.eecs.umich.edu/~msmelyan/papers/physsim_onmanycore_itj.pdf
[Box2D - SIMD Matters]: https://box2d.org/posts/2024/08/simd-matters/
[Erin Catto]: https://github.com/erincatto

## 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:

![Large Pyramid 2D benchmark](https://github.com/user-attachments/assets/6cd1b6c0-a59e-4b0a-951d-f0e677a1d330)

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:

![Many Pyramids 2D benchmark](https://github.com/user-attachments/assets/8f315fac-ef8b-421a-8a75-db6c593feb82)

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:

![Before](https://github.com/user-attachments/assets/09d9019d-3066-458f-a1cc-16bd4918b7fa)

And now with graph coloring, on this branch:

![After](https://github.com/user-attachments/assets/de00d88c-e868-40aa-805f-49e7312beee0)

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 `ContactEdge`s in the graph, containing only edge connectivity information with stable `ContactId`s.
- The `ContactGraph` now stores `ContactPair`s 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.