ODE Simulation Stability         Changing the way ODE is called, changes the behaviour of your simulation
Image

Simulation Stability


Changing the way ODE is called, changes the behaviour of your simulation. As far as possible, I’d recommend getting your physics ‘behaving well’ before you spend too long tweaking your objects physical properties.

I started with the Demo_BSPCollision sample application, and made minor changes like increasing the height of the box stack, and changing to the ‘shelf’ to a finite plane. I discovered that the simulation suffered from jittering, meaning the boxes would never fully settle, and occasionally would jump around unpredictably.

Looking at the forums, this is not uncommon. The regular solutions are to shorten the simulation step time, and to increase the amount of collision two objects can have with each other.

Note: I've duplicated and renamed the ref app classes in building my program. I have attempted to rename things back to the original naming for this article, but I may have missed a few things. The code may not compile straight away. Feel free to correct it further.

Debugging and step size

My first problem was that my old super-slow machine would lockup in debug. Each physics update would take longer to calculate than the time it was simulating, which put it in a viscious circle of having to simulate larger and larger time slices, and therefore taking longer and longer to execute. After just a few loops, my machine was pretty much locked up. Another problem was that if I stopped in the debugger for a while, when I continued executing the simulation would then simulate all the time that I had spent in the debugger.

These can both be fixed by clamping the maximum timestep for a single frame. There are likely nicer ways to do this, but here’s what I did:

<nowiki>void World::simulationStep(Real timeElapsed)
{
        /* Hmm, gives somewhat jerky results*/
        static Real leftOverTime = 0.0f;

    Real time = timeElapsed + leftOverTime;    

    // HACK – Clamp maximum step size
    if(time > 0.2f)
        time = 0.2f;</nowiki>

Increasing Simultaneous Collisions

It’s pretty easy to change ‘ApplicationObject::testCollide’ to handle N collisions instead of just one. But to save you the time, here’s how mine now looks:

<nowiki>bool ApplicationObject::testCollide(ApplicationObject* otherObj)
{
    const int MAX_COLLISIONS = 4;
    int collisionCap = MAX_COLLISIONS;    // For comparing behaviours
    bool collided = false;
    dContactGeom contactGeomList[MAX_COLLISIONS];
    dContactGeom *pContactGeom;
    dGeom *o1, *o2;
    CollisionProxyList::const_iterator proxy1, proxy2, proxy1end, proxy2end;
    proxy1end = mCollisionProxies.end();
    proxy2end = otherObj->mCollisionProxies.end();

    CollisionInfo collInfo;

    for (proxy1 = mCollisionProxies.begin(); proxy1 != proxy1end; ++proxy1)
    {
        for (proxy2 = otherObj->mCollisionProxies.begin(); proxy2 != proxy2end; ++proxy2)
        {
            o1 = *proxy1;
            o2 = *proxy2;
            int numc = dCollide(o1->id(), o2->id(), MAX_COLLISIONS, (dContactGeom*)&contactGeomList, sizeof(dContactGeom));

            pContactGeom = (dContactGeom*)&contactGeomList;
            if(numc > collisionCap)
                numc = collisionCap;
            while (numc)
            {
                // Create contact joints if either object is dynamics simulated
                // If one is not, then sim will not affect it anyway, it will be fixed
                // However if one is enabled, we need the contact joint
                if (this->isDynamicsEnabled() || otherObj->isDynamicsEnabled())
                {
                    // We use the most agressive parameters from both objects for the contact
                    dContact contact;
                    Real bounce, velThresh, softness;
                    // Use the highest coeff of restitution from both objects
                    bounce = std::max(this->getBounceRestitutionValue(), 
                            otherObj->getBounceRestitutionValue());
                    // Use the lowest velocity threshold from both objects
                    velThresh = std::min(this->getBounceVelocityThreshold(),
                            otherObj->getBounceVelocityThreshold());
                    // Set flags
                         contact.surface.mode = dContactBounce | dContactApprox1;
                    contact.surface.bounce = bounce;
                    contact.surface.bounce_vel = velThresh;

                    softness = this->getSoftness() + otherObj->getSoftness();
                    if (softness > 0)
                    {
                        contact.surface.mode |= dContactSoftCFM;
                        contact.surface.soft_cfm = softness;
                    }
                        
                    // Set friction to min of 2 objects
                    // Note that ODE dInfinity == Math::POS_INFINITY
                    contact.surface.mu = std::min(this->getFriction(), otherObj->getFriction());
                    contact.surface.mu2 = 0;
                    contact.geom = *pContactGeom;
                    dContactJoint contactJoint(
                            World::getSingleton().getOdeWorld()->id(), 
                            World::getSingleton().getOdeContactJointGroup()->id(), 
                            &contact);

                    // Get ODE bodies
                    // May be null, if so use 0 (immovable) body ids
                    dBody *b1, *b2;
                    dBodyID bid1, bid2;
                    bid1 = bid2 = 0;
                    b1 = this->getOdeBody();
                    b2 = otherObj->getOdeBody();
                    if (b1) bid1 = b1->id();
                    if (b2) bid2 = b2->id();
                    contactJoint.attach(bid1, bid2);
                }

                // Tell both objects about the collision
                collInfo.position.x = pContactGeom->pos[0];
                collInfo.position.y = pContactGeom->pos[1];
                collInfo.position.z = pContactGeom->pos[2];
                collInfo.normal.x = pContactGeom->normal[0];
                collInfo.normal.y = pContactGeom->normal[1];
                collInfo.normal.z = pContactGeom->normal[2];
                collInfo.penetrationDepth = pContactGeom->depth;
                this->_notifyCollided(otherObj, collInfo);
                otherObj->_notifyCollided(this, collInfo);

                // set return 
                collided = true;

                pContactGeom++;
                numc--;
            }
        }
    }
    return collided;
}</nowiki>


When I first increased the number of collisions I was asking for, I started to get a stack overflow with a large box stack. Looking at the ODE help it says the following:

ODE with dWorldStep requires stack space roughly on the order of O(n)+O(m2), where n is the number of bodies and m is the sum of all the joint constraint dimensions. If m is large, this can be a lot of space!
Unix-like operating systems typically allocate stack space as it is needed, with an upper limit that might be in the hundreds of Mb. Windows compilers normally allocate a much smaller stack. If you experience crashes when running large systems, try increasing the stack size. For example, the MS VC++ command line compiler accepts the /Stack
num flag to set the upper limit.
Another option is to switch to dWorldQuickStep.


Information from here:

I first tried putting my stack size up to over 100 meg (from the .net default of 1 meg) but the problem persisted for me. So I instead decided to change to dWorldQuickStep.

I’ve read on the Ogre forum that Ogre uses the ODE ‘dWorldQuickStep’ stepping method, instead of the more complex ‘dWorldStep’ method. Stepping through the code however the reference app framework uses dWorldStep (in Ogre V 1.0.1 at least).

If you wish to switch to the dWorldQuickStep, and are building on the ref app code, you can change OgreRefAppWorld.cpp (or your copy of it) as follows. Replace:

mOdeWorld->step(dReal(timeElapsed));

with:

dWorldQuickStep( mOdeWorld->id(), dReal(timeElapsed) );

You could instead just make the change inside the step function, but for my purposes it was ‘cleaner’ to do it outside.

After changing to dWorldQuickStep and increasing the number of simultaneous collisions from one to four, my simulation looked noticeably more stable, but still contained some small jittering in debug and very small jittering in release. The ‘very small’ jittering in release isn’t a big problem to look at, but ideally it should become static and stop performing physics calculations. The fact that release is better than debug is due to debug taking larger time steps.

Reducing Step Size

My simulation is running with reasonable stability right now, so I’m going to play with other things for a while. I’ve not played with step size much yet, but you might want to call this:

mWorld->setSimulationStepSize(0.02f);

And inspect the simple inner workings of this function:

void World::simulationStep(Real timeElapsed)

ERP - Error Reduction Parameter



ODE has a global variable called the Error Reduction Parameter (ERP). Occasionaly the physics objects will get into a situation where they're not quite following the rules - joints don't line up correctly, or rigid objects might be slightly inside one another. The ERP parameter specifies how quickly ODE is to fix these errors. Solving them quickly can result in jittery behaviour, solving them slowly smoothes things out, but may mean the user can see the errors gradually getting fixed. To modify the ERP you can either use the following function from odecpp.h or directly call the function that it wraps. For more information on ERP see the ODE documentation.

class dWorld {
 ...
  void setERP (dReal erp)
    { dWorldSetERP(_id, erp); }

Postscript

After writing the page above I still had some instability. A bigger problem for me was that I realized that (to my understanding) ODE was non-deterministic. This means that if you try and run the same simulation twice, you get slightly different different results. If you were, say, just throwing some pieces of crate around for an effect, this would not be an issue. My needs however required determinism.

After hearing positive praise about it, I switched over to Newton physics. I found it to be very stable. The code interface is a little different, but good. Newton claims to be deterministic, and for the most part it is very close to being, though using test cases with many objects can show it to be non deterministic. It seems good enough for my needs. The porting process took me a little while, but the sample applications helped.

Overall, my personal preference is now Newton. I apologize if these comments are not entirely correct, it may be the case that I've not got the best out of ODE. Also, if I hadn't had the determinism issue, ODE would have been fine...