Solvers in bullet

[KenB]

Hi,
I notice there are several different solvers in Bullet e.g. borrowed from ODE and also a solver based on Catto's work. What is the overall status of these solvers in comparison and which one is generally used?
Are they all "working and stable" or is there "work in progress" in any of them?

[Erwin Coumans]

The default constaint solver in Bullet, btSequentialImpulseConstraintSolver, was developed in 2003. A few years ahead of Erin Catto's presentations, although he pointed out some bug related to clamping accumulated impulses.

What is the overall status of these solvers in comparison and which one is generally used?

They are essentially the same methods, sequential impulse boils down to projected gauss seidel (PGS). btSequentialImpulseConstraintSolver is default, ODE quickstep solver (PGS) is just optional in the Extras folder. This reminds me of a 1999 paper by Moreau, that Kenny Erleben brought to my attention recently.
Numerical aspects of the sweeping process refers to Gauss Seidel in the rigid body/contact context:

9.2: Here is an iteration technique a la Gauss-Seidel which consists of treating a succession of single-contact problems.

btSequentialImpulseConstraint is the default choice. It has been further optimized and parallelized for Playstation 3 Cell SPUs and XBox 360 Xenon/multi-core. Our latest SPU optimizations exceed Ageia PhysX performance on Cell. The quickstep solver is just an optional extra, albeit currently faster/more cache-friendly on PC.

Are they all "working and stable" or is there "work in progress" in any of them?

They should be working and reasonably stable, but they suffer common issue with iterative solvers: difficulty of handling large mass ratios, when stacking a very heavy object on top of a light object. We could add some solver quality improvements, which would make it better then ODE's quickstep, see Erin's posting for some ideas.

The development is mainly driven by feedback, so if you have any suggestions please let it know,
Thanks,
Erwin

[Dirk Gregorius]

Numerical aspects of the sweeping process refers to Gauss Seidel in the rigid body/contact context:

Cool, thanks for the paper. This is one of the few I was missing

[KenB]

Thanks for the reply Erwin. Yes, I agree that they are all similar in the sense that they use projected Gauss-Seidel but the specific formulations aren't available outside the code it seems.

The specifics in how friction is implemented, how clamping is done, the balance of impact velocity vs penetration dependent impulses and how this all combines into under- or over-relaxation, can make a big difference at the end of the day, and right now these things are deep secrets hidden in the code

I just noticed that the contact generation in Bullet delivers one contact at a time and therefore e.g. a
box bouncing flat on a plane gets a huge symmetry breakage. Boxes typically don't fall flat (compare a dice) so this isn't much of a problem, but it gives some nifty problems when simulating e.g. a sledgehammer that wants to bounce off its constraint due to the symmetry breakage. The combination of low order integration, large timestep, and few iterations - with an assymmetric contact didn't quite give the result I was looking for...
Is this a fundamental thing with contact generation in Bullet?

[Dirk Gregorius]

The specifics in how friction is implemented, how clamping is done, the balance of impact velocity vs penetration dependent impulses and how this all combines into under- or over-relaxation, can make a big difference at the end of the day, and right now these things are deep secrets hidden in the code

I think Bullet (and Erwin might correct me when I am wrong) does the following:

1) Friction is implemented as two linear constraints in the contact plane at each contact point and clamped against the associated normal impulse. In order to improve friction the first friction direction is aligned with the relative velocity in the tangent plane.

2) I think there is no special code to balance the impact velocity and the penetration dependent impulse. But I am not really sure what you mean. Are you only refering to contacts or arbitrary constraints? In formulas Bullet does: lambda = -( tau * C / dt - J * v ) / (J * W * JT ) and P = JT * lambda ( tau usually 0.1 - 0.2 ).

3) Not sure, but I believe Bullet doesn't use under- or overrelaxation. The experience shows that overrelaxation adds energy if you exit the iteration process early like we do for games. It might be good if you can wait for convergence. I never tested this.


Maybe you can share some of your experiences . Especially I would be interested if you suggest different friction models.


HTH,
-Dirk

[Erwin Coumans]

Thanks for the reply Erwin. Yes, I agree that they are all similar in the sense that they use projected Gauss-Seidel but the specific formulations aren't available outside the code it seems.

The specifics in how friction is implemented, how clamping is done, the balance of impact velocity vs penetration dependent impulses and how this all combines into under- or over-relaxation, can make a big difference at the end of the day, and right now these things are deep secrets hidden in the code

I'm in the process of writing a book that should cover every detail, with tentative title 'Advanced Game Physics' for the Morgan Kaufmann 3D graphics series. Right now is most busy time of the year, preparing for Game Developers Conference next week, so please let me get back to this after GDC.

I'll try to get back to those code snippets with line numbers from btSequentialImpulseConstraintSolver.cpp.

Don't use penetration depth recovery, if the restitution already recovers the same amount:

Code: Select all

							btScalar penVel = -solverConstraint.m_penetration/infoGlobal.m_timeStep;
  665 							solverConstraint.m_penetration *= -(infoGlobal.m_erp/infoGlobal.m_timeStep);
  666 
  667 							if (solverConstraint.m_restitution > penVel)
  668 							{
  669 								solverConstraint.m_penetration = btScalar(0.);
  670 							}

Friction model setup:

Code: Select all

 690 							btVector3 frictionDir1 = vel - cp.m_normalWorldOnB * rel_vel;
  691 							btScalar lat_rel_vel = frictionDir1.length2();
  692 							if (lat_rel_vel > SIMD_EPSILON)//0.0f)
  693 							{
  694 								frictionDir1 /= btSqrt(lat_rel_vel);
  695 								addFrictionConstraint(frictionDir1,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,colObj0,colObj1, relaxation);
  696 								btVector3 frictionDir2 = frictionDir1.cross(cp.m_normalWorldOnB);
  697 								frictionDir2.normalize();//??
  698 								addFrictionConstraint(frictionDir2,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,colObj0,colObj1, relaxation);
  699 							} else
  700 							{
  701 								//re-calculate friction direction every frame, todo: check if this is really needed
  702 								btVector3	frictionDir1,frictionDir2;
  703 								btPlaneSpace1(cp.m_normalWorldOnB,frictionDir1,frictionDir2);
  704 								addFrictionConstraint(frictionDir1,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,colObj0,colObj1, relaxation);
  705 								addFrictionConstraint(frictionDir2,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,colObj0,colObj1, relaxation);
  706 							}
  707 
  708 						}

I just noticed that the contact generation in Bullet delivers one contact at a time
[...]
Is this a fundamental thing with contact generation in Bullet?

We should provide two improved methods, so users can choose higher quality contact generation.
The current default of adding one point at a time is just one option that seems to work fine in many cases. But there are exceptions as you already pointed out.

See also the discussion related to this here:
http://www.bulletphysics.com/Bullet/php ... f=4&t=1857

I added an issue in Bullet.GoogleCode tracker, and follow up on this later:
See http://code.google.com/p/bullet/issues/detail?id=20

Thanks,
Erwin