Rework contact solver and collision detection, implement speculative collision

[Jondolf]

Objective

Note: This is for the avian branch, which will later be pushed to main.

Bevy XPBD currently uses Extended Position-Based Dynamics for its solver. While it has worked quite well, it does have several issues:

• Solving deep overlap is very energetic and explosive.
• Position-based constraints often have high-frequency oscillation where bodies are always slightly moving back-and-forth. XPBD rarely (if ever) reaches a truly stable and relaxed state without sleeping.
• A ton of substeps are needed for minimizing jitter and making systems of bodies stable.
• Friction is typically not quite as robust as it is with impulse-based methods. There are several ways of handling it, some better than others, but it can still have some issues.
• Many features can be more challenging to implement and expose APIs for, and there are fewer learning resources and references for the usage of XPBD in game physics.

Additionally, XPBD even has some potential legal ambiguities, see #346.

Aside from XPBD itself, there are also several other issues with the current collisions and rigid body dynamics:

• Narrow phase collision detection is run at every substep instead of just once per frame like most physics engines. This can be very bad for performance.
• There is no support for any kind of Continuous Collision Detection, so tunneling can be a big issue.
• It is not possible to configure the behavior of collisions aside from the number of substeps.
• Collision events are sent when not actually in contact #224
• Sleeping in 0.4 just doesn't work.

These are big issues, and I don't currently see XPBD as the best option for the contact solver going forward. Extensive rewrites and fixes are needed.

Solution

This PR consists of three major parts:

• Rewrite the contact solver to use an impulse-based TGS Soft solver instead of XPBD.
• Rework the narrow phase and other collision detection logic.
• Implement speculative collision, a form of Continuous Collision Detection.

They all relate to each other, and are covered in detail below.

Substepped Impulse-Based Solver With Soft Constraints

The contact solver has been rewritten to use TGS Soft, a substepped impulse-based solver using soft constraints. The choice was largely motivated by the wonderful Erin Catto's Solver2D experiments, where TGS Soft was deemed as quite a clear winner.

Box2D V3 was used as the primary inspiration for the core implementation of the new solver. Engines such as Rapier and Bepu also use a very similar approach.

Terminology

• Projected Gauss-Seidel (PGS): The classic iterative approach to solving constraints (contacts and joints) using the Gauss-Seidel numerical method. Reframed by Erin Catto as Sequential Impulses.
• Temporal Gauss-Seidel (TGS): Same as PGS, but prefers substepping over iteration, running more simulation steps with smaller time steps rather than simply iteratively solving constraints. Substeps tend to be more effective than iterations, as shown in the Small Steps in Physics Simulation paper by Macklin et al.
• Baumgarte stabilization: When solving contact constraints using impulses, boost the impulses using a bias to account for overlap and actually push the bodies apart.
• Soft constraints: Similar to Baumgarte stabilization, but more stable and controlled. Based on the harmonic oscillator, soft constraints dampen constraint responses, and can be tuned intuitively with a frequency and damping ratio.
• Warm starting: Store the constraint impulses from the previous frame and initialize the solver by applying them at the start of the current frame. This helps the solver converge on the solution faster, and is especially helpful when objects are coming to rest.
• Relaxation: Baumgarte stabilization and soft constraints can add unwanted energy. Relaxation helps remove it by solving constraints a second time, but without a bias.

Please refer to the Solver2D post for a more complete overview of what TGS and soft constraints are, how they work, and how they relate to other approaches.

Solver Overview

As stated earlier, contacts use TGS Soft. However, joints still currently use XPBD, so the new solver is actually a kind of hybrid solver. I do plan on transitioning joints away from XPBD in the future though.

Below is a high-level overview of the new structure of the solver.

1. Broad phase collision detection collects potential collision pairs into BroadCollisionPairs.
2. Narrow phase collision detection computes contacts for the pairs and adds them to Collisions.
3. A ContactConstraint is generated for each contact manifold, and added to ContactConstraints.
4. Substepping loop, running SubstepCount times.
   1. Integrate velocities, applying gravity and external forces.
   2. Warm start the solver.
   3. Solve velocity constraints with bias (soft constraints).
   4. Integrate positions, moving bodies based on their velocities.
   5. Relax velocities by solving constraints again, but without bias.
   6. Solve XPBD constraints (joints) and perform XPBD velocity projection.
5. Apply restitution as a post-phase.
6. Finalize positions by applying AccumulatedTranslation.
7. Store contact impulses for next frame's warm starting.

Refer to the code for implementation details. The contact logic and constraints are quite heavily commented and should hopefully be relatively easy to follow.

New Solver Results

Collisions have significantly less drift than before, and performance is much better. Below is a pyramid with a base of 64 boxes, simulated with 4 substeps, with sleeping disabled.

Old: XPBD has substantial drift, and the pyramid quickly collapses in on itself. This even happens with a much larger number of substeps, although to a slightly lesser degree. Performance is very poor, even with just 4 substeps.

2024-06-24.22-01-47.mp4

New: With TGS Soft, the pyramid stays stable. There is a very small amount of drift over a long period of time, but even that can be mitigated by configuring the contact behavior through the SolverConfig and/or by adding more substeps.

2024-06-24.21-58-54.mp4

Impulses even out and stabilize much better with TGS Soft, even with deep overlap. In overlapping cases, the old XPBD implementation was significantly more explosive. Below is an example where colliders are dynamically enlarged to fill up a container.

Old: With XPBD, overlap causes explosions even before the shapes fill the container. Once they do fill the container, they jitter very strongly and explode through the walls. This is typically very undesirable for games.

2024-06-24.22-20-27.mp4

New: Overlap is solved perfectly smoothly with no explosions. The contact impulses even out without jitter. With no space to go, the shapes prioritize stability over perfectly hard contacts that would cause explosiveness.

2024-06-24.22-22-06.mp4

An important thing to note is that the new solver uses several tolerances and thresholds that are length-based. The old solver also had some tolerances, but now they are much more important for stability.

Without any tuning, a 2D game using pixels as length units might have collision issues, because the tolerances would be wrong for that application. For example, below is a scene with stacks of balls that have a radius of 5.0, with no tuning whatsoever:

The contacts are too soft. To fix this, there is a new PhysicsLengthUnit resource, which can be thought of as a kind of pixels-per-meter conversion factor. It is only for scaling the internal tolerances (and debug rendering gizmos!) however, and doesn't scale objects or velocities in any way.

PhysicsLengthUnit can be easily set when adding PhysicsPlugins for an app:

fn main() {
    App::new()
        .add_plugins((
            DefaultPlugins,
            // A 2D game with 20 pixels per meter
            PhysicsPlugins::default().with_length_unit(20.0),
        ))
        .run();
}

And with that, we get the stability and behavior we expect:

Collision Detection Refactor

The narrow phase and the PostProcessCollisions schedule are now run in the PhysicsSchedule instead of the SubstepSchedule. This greatly improves performance, and also makes collision events and reacting to collisions less awkward and footgunny: previously, the contact data was always from the last substep, but at that point, the contact is often already mostly resolved.

Moving narrow phase collision detection out of the substepping loop is possible by simply storing contact data in local space and computing the updated separation distance at each substep using the current poses of the bodies. This way, the contact data is (approximately) accurate even with the bodies moving relative to each other within substeps.

A few other improvements have been made as well:

• The narrow phase logic has been extracted into a NarrowPhase system parameter.
• The broad phase now outputs intersection pairs with the entities in ascending order instead of based on the minimum X extent. This fixes issues where bodies moving past each other on the X axis counts as a separate collision, causing problems with collision events.
• Collision events are sent when not actually in contact #224 is fixed.
• The warnings logged for overlapping bodies have less false alarms.

Speculative Collision

Tunneling is a phenomenon where fast-moving small objects can pass through thin geometry such as walls due to collision detection being run at discrete time steps:

Moving the narrow phase out of the substepping loop has the unfortunate consequence that it increases the risk of tunneling as collisions are not computed as frequently. A solution for this is needed before we commit to the narrow phase change.

One of the primary solutions to tunneling is Continuous Collision Detection. There are two common forms:

1. Sweep-based CCD: Each collider with CCD enabled is swept from the current position to the predicted one (or from the previous position to the current one). If a hit is detected, the bodies are moved to the time of impact. Contact resolution can be left to the next frame, or to avoid "time loss", performed immediately with a substepping scheme that can also consider secondary collisions.
2. Speculative CCD: Contact points are predicted by the narrow phase for fast-moving objects before they actually touch. Speculative collision response is handled by pushing back the part of velocity that would cause penetration.

Sweep-based CCD can be much more expensive and more complex, especially with substepping and non-linear sweeps. Speculative collision on the other hand is very cheap while still being quite robust, although in rare cases it can miss collisions or cause ghost collisions.

The new solver implements speculative collision, which is enabled for all bodies by default. In my experiments, this is efficient, improves stability, and eliminates almost all tunneling except when caused by contact softness. This approach seems to also be taken by both Bepu and Box2D.

Below is a high-level overview of how speculative collision is currently implemented.

• The speculative margin is the maximum distance at which a collision pair generates speculative contacts. This is unbounded for every rigid body by default, which eliminates almost all tunneling.
• Each AABB is expanded in the movement direction based on the body's velocity, clamped by the speculative margin if it is bounded. (Box2D might have a way to avoid this expansion, but I haven't looked into it yet)
• The effective speculative margin is the actual margin used for contact computation, and it is clamped based on the velocities. This is used as the maximum separation distance for contact / closest point computation.
• When actually solving the contact (normal part, no friction), use the softness parameters only if the contact is penetrating. Otherwise it is speculative, and we bias the impulse to cancel out the velocity that would cause penetration.
• Apply restitution in a separate phase after the substepping loop.

From a user's point of view, this happens completely in the background. However, if desired, the speculative margin can be configured for an entity using the SpeculativeMargin component, or even globally using default_speculative_margin in NarrowPhaseConfig.

Below is an example of balls being shot at high speeds at thin walls and spinning objects, with just a single substep to eliminate the effect of substepping. The simulation is paused and stepped manually a few times to closer see the behavior.

Old: Almost all of the balls simply pass through the thin geometry. With more substeps, this tunneling could been reduced, but never removed completely. Increasing the substep count would also hurt performance.

2024-06-25.17-34-10.mp4

New: The balls never tunnel straight through the walls, and often even hit the spinning objects. You can still see many balls phasing through the walls, but this is mostly due to contact softness, in cases where a body hits another body hard enough to force it through the wall.

2024-06-25.17-30-44.mp4

Of course, bodies getting pushed through the ground by other bodies is still an issue. This could be reduced in the future by solving contacts against static objects after dynamic-dynamic contacts, giving them higher priority and stiffness. This is follow-up material however, and cases where this is an issue should be quite rare in games.

Sweep-based CCD will also be added as an option in a follow-up. I already have an implementation ready locally.

Performance Results

Comparing this branch to the avian branch (which still has the old solver), there is a roughly 4-5x performance improvement for collision-heavy scenes, with the difference growing with the number of collisions and the number of substeps.

The benchmarks below use 8 substeps.

With just a single substep, the difference is smaller, but this branch is still faster, up to 2x.

Other Changes

This is such a big rework that it is unfortunately almost impossible to cover every change. Below are some noteworthy changes however.

Sleeping Rework

• Sleeping in 0.4 was broken, and bodies never fell asleep properly. This has been fixed.
• Bodies now store a PhysicsChangeTicks component. This is used to detect what component changes are made by the user outside of physics schedules, and what changes are made by the physics engine. This way, we can ignore changes made by the physics engine and have more control over when bodies are woken up.
• For now, bodies can only sleep when they are not in contact with other dynamic bodies. This is because the current per-body sleeping approach is quite buggy and has stability issues. In the future, I plan on implementing simulation islands and doing per-island sleeping like most other physics engines.

Integrator Rework

• Velocity integration and position integration now run in separate systems and in different parts of the schedule: IntegrationSet::Velocity and IntegrationSet::Position. This is needed for the new solver.
• The semi-implicit Euler integration scheme now has its own module with proper docs and tests.
• Integration uses par_iter_mut.

Scheduling Changes

• PhysicsStepSet now has First and Last variants, so users can easily schedule systems before or after physics in the PhysicsSchedule.
• The narrow phase and the PostProcessCollisions schedule are now run in PhysicsStepSet::NarrowPhase instead of SubstepSet::NarrowPhase.
• Integration is now run in IntegrationSet::Velocity and IntegrationSet::Position instead of SubstepSet::Integrate.
• SubstepSet has been removed.
  • The solver runs in PhysicsStepSet::Solver.
  • The solver's system sets are in SolverSet.
  • Substepping is performed in SolverSet::Substep.
  • The substepping loop's system sets are in SubstepSolverSet.

New Configuration Options

• NarrowPhaseConfig has new default_speculative_margin and contact_tolerance properties.
• The new SolverConfig resource can be used for tuning collisions.
• The new PhysicsLengthUnit resource can be used as a scaling factor for internal length-based tolrances.
  • Existing length-based tolerances and thresholds are now scaled by this, like the linear property of SleepThreshold.
  • Debug rendering is also scaled by the length unit.

Examples

• XpbdExamplePlugin has been renamed to ExampleCommonPlugin, and it no longer adds PhysicsPlugins automatically.
• 2D examples have PhysicsLengthUnit configured.
• The one_way_platform example's logic has been modified to account for PostProcessCollisions no longer running in the SubstepSchedule.
• The collision logic for the kinematic character controller examples has been rewritten to run in the PhysicsSchedule without jitter or stability issues.

Miscallaneous

• ColliderAabb now has the grow and shrink methods.
• ContactData now stores feature IDs for contact matching, which is needed for warm starting.
• ContactData property index has been removed.
• 3D tangent impulses are now 2D vectors instead of scalar values.

---

Migration Guide

New Contact Solver

The contact solver has been rewritten. In practice, this has the following effects:

• Collisions should be much more stable.
• Resolving overlap is no longer nearly as explosive.
• Less substeps are generally needed for stability.
• Tunneling is much more rare.
• Performance is better.

However:

• Contacts may even be too soft by default for some applications. This can be tuned with the SolverConfig.
• Static friction is currently not considered separately from dynamic friction. This may be fixed in the future.
• Restitution might not be quite as perfect in some instances (this is a tradeoff for speculative collision to avoid tunneling).
• 2D applications may need to configure the PhysicsLengthUnit to get the best stability and behavior.

The PhysicsLengthUnit can be thought of a pixels-per-meter scaling factor for the engine's internal length-based tolerances and thresholds, such as the maximum speed at which overlap is resolved, or the speed threshold for allowing bodies to sleep. It does not scale actual physics objects or their velocities.

To configure the PhysicsLengthUnit, you can insert it as a resource, or simply set it while adding PhysicsPlugins:

fn main() {
    App::new()
        .add_plugins((
            DefaultPlugins,
            // A 2D game with 20 pixels per meter
            PhysicsPlugins::default().with_length_unit(20.0),
        ))
        .run();
}

Custom Constraints

Custom constraints using XPBD are currently still possible. However, the traits and systems are now located in the dynamics::solver::xpbd module, and users should run solve_constraints in SubstepSolverSet::SolveUserConstraints instead of SubstepSet::SolveUserConstraints.

Scheduling Changes

Several scheduling internals have been changed. For example:

• The narrow phase and PostProcessCollisions schedule are now run in PhysicsStepSet::NarrowPhase instead of SubstepSet::NarrowPhase.
• Integration is now run in IntegrationSet::Velocity and IntegrationSet::Position instead of SubstepSet::Integrate.
• SubstepSet has been removed.
  • The solver runs in PhysicsStepSet::Solver.
  • The solver's system sets are in SolverSet.
  • Substepping is performed in SolverSet::Substep.
  • The substepping loop's system sets are in SubstepSolverSet.

Systems running in PostProcessCollisions may need to be modified to account for it being moved outside of the substepping loop.

[janhohenheim]

I think a potential new user will be a bit confused by this line.
This goes for the other instances of with_length_unit in the other examples as well. Is there a reason why the examples don't use the default unit length?

[Jondolf]

These examples have a 2D camera with a default projection, so the objects use pixels for their dimensions and are quite large in terms of physics units if the length unit was the default 1.0. This leads to the issues with overly soft contacts and incorrectly tuned tolerances, which can be seen at the end of the "New Solver Results" section (the two screenshots).

Note that bevy_rapier also has this, but in 2D it calls the method pixels_per_meter, even though internally the property is still called length_unit. And it also has a with_length_unit method, so it just has two differently named methods that essentially configure the same thing.

A quick survey on what names different engines use for this property:

• Rapier: length_unit (or in 2D, also the pixels_per_meter method)
• Box2D: b2_lengthUnitsPerMeter
• PhysX: PxTolerancesScale.length

I agree though that the method usage should probably be commented in the examples, or the method name itself could be more self-explanatory. It's just hard to come up with an approachable name that conveys the purpose, but does not imply that it actually scales velocities or forces in any way.

For example, if there was a physics_scale property and it was set to 10.0 (let's imagine this is 2D), a user could reasonably assume that Collider::rectangle(1.0, 1.0) would get scaled to a 10x10 pixel collider. Or you could assume that you could use LinearVelocity::X, and have the body move at the speed of 10 pixels per second. But this is not what length_unit does, as it just scales the internal length-based tolerances and thresholds to match the expected scale. The nuance is difficult to convey through the name of the property.

For these examples though, another option would be to just scale the camera projection so that physics objects don't have to be so large in terms of physics units and you don't need a custom length unit. This would technically be more ideal for physics in other ways too, although I don't know what the game dev best practices are for 2D games when it comes to using pixels vs. other units for scenes.

(sorry for the long response lol)

[Jondolf]

For now, I added a comment briefly explaining the length unit to each example that sets it. I also improved the length unit docs overall

[janhohenheim]

Too bad minmax is nightly only :/

[Jondolf]

Addressed most of the comments. A few other changes:

• IntegratorPlugin initializes its system sets in the schedule given to it instead of only configuring them in SolverPlugin. This way, it would work on its own even without the solver.
• Added a match_contacts property to NarrowPhaseConfig. Previously it was decided by SolverConfig::warm_start_coefficient, but that's a solver implementation detail.
• NarrowPhasePlugin has a generate_constraints property. It's unnecessary to generate constraints if you only use the narrow phase for collision detection without the solver.
• NarrowPhasePlugin can be initialized with any schedule, not just PhysicsSchedule.