When I implement a 3d physics engine, stacking box fall easily. ThenI added warm start, it's more stable now. But it still slightly slide to one side when 10 or more boxes stack.

I don't know if there are bugs in my physics engine, or just precision or parameter fine-tuning problem. Can somebody help? thanks sincerely.

Some things that may help:

• Randomize the order of pairwise constraint solving on each solver iteration. This prevents the same constraints from always dominating the results since they are applied last, and helps to propagate information from the top to bottom of stack.
• If you aren't already, a fixed time step also can improve stability.
• Another thing I have done is to first partition the graph of all constraints into disjoint parts before solving. The graph is split at static objects into constraint “islands” (actually built in a bottom-up clustering way). Then solve each island independently. With this splitting of the graph, you can do a different number of solver iterations in each island. In my implementation, I stop solver iterations in an island if the largest impulse applied is below a threshold. This means most small islands of a few constraints are solved quickly in 1 or 2 iterations, but more complex stacks/chains automatically get more iterations for stability.
• Make sure your contact generation is accurate (including positions, normals, penetration depth). You should be getting at least 3-4 contacts per pair of touching boxes. It can also help to treat each contact as a pair of collision points (one on body A and one on body B), rather than at a single point.

Aressera said:

Some things that may help:

• Randomize the order of pairwise constraint solving on each solver iteration. This prevents the same constraints from always dominating the results since they are applied last, and helps to propagate information from the top to bottom of stack.
• If you aren't already, a fixed time step also can improve stability.
• Another thing I have done is to first partition the graph of all constraints into disjoint parts before solving. The graph is split at static objects into constraint “islands” (actually built in a bottom-up clustering way). Then solve each island independently. With this splitting of the graph, you can do a different number of solver iterations in each island. In my implementation, I stop solver iterations in an island if the largest impulse applied is below a threshold. This means most small islands of a few constraints are solved quickly in 1 or 2 iterations, but more complex stacks/chains automatically get more iterations for stability.
• Make sure your contact generation is accurate (including positions, normals, penetration depth). You should be getting at least 3-4 contacts per pair of touching boxes. It can also help to treat each contact as a pair of collision points (one on body A and one on body B), rather than at a single point.

Thanks for your nice advice. @Aressera

1. If I randomize the order of pairwise constraint, will it introducing non-determinism? I'll try this.
2. a fixed time step is already used at 60Hz
3. I do not build split island to solve constraints. As you mentioned, more complex stacks/chains automatically get more iterations for stability. In my implementation, I iterate 30 times for solving velocity, 10 times for solving position. I compared this with Jolt physics, Jolt performs better than me at same iterate times.
4. I used contact manifold. when two boxes contact, there are usually 4 points to build contact constraint. But I noticed, after several minutes passed, box 1 and 2 (the lowest 2 boxes) have only 3 contact point. one point in manifold's distance to another box's face is larger than 0.02m. So I filtered this point from manifold.