Open Dynamics Engine

The Open Dynamcis Engine (herein referred to as 'ODE') is a physics engine intended for computer games scientific research. From the site:

ODE is an open source, high performance library for simulating rigid body dynamics. It is fully featured, stable, mature and platform independent with an easy to use C/C++ API. It has advanced joint types and integrated collision detection with friction. ODE is useful for simulating vehicles, objects in virtual reality environments and virtual creatures. It is currently used in many computer games, 3D authoring tools and simulation tools.

While initial tests proved promising, the engine was ultimately found lacking, as described below. The Bullet Engine may serve as a better alternative, though testing has not been performed to indicate this one way or the other.

Basic ODE Implementation

In order to test ODE for our purposes, we performed both sphere and polyhedral simulations using the software and compared the sphere-simulation output with known correct output in order to verify the accuracy of the engine.

Sphere Simulations

The original ODE testing was done with spheres using the equations given in Local Simulations of Planetary Rings (Wisdom and Tremaine, The Astronomical Journal, Vol. 95, Num. 3, March 1988). The results seemed promising, and guiding center motion was clearly preserved (see side image). The code for this simulation is relatively straightforward, and is fairly easy to work with.

The first piece of user-defiend code is the collisionCheck function which takes two geometry objects and a 'data' argument. Below I get the geometry body died to the body id passed in and create a 'contact' object which represents a potential collision in ODE. I set the restitution of the surface (called bounce in the code) and define its friction (mu in the code) as infinity for testing purposes. Note that in order to process a collision, ODE does the following: first, check for collisions, second, create joints for any collisions, and third, add these joints to a group of joints to be processed that represents all collisions.

// this is called by dSpaceCollide when two objects in space are potentially colliding
static void collisionCheck(void *data, dGeomID o1, dGeomID o2) {
   dBodyID b1 = dGeomGetBody(o1);
   dBodyID b2 = dGeomGetBody(o2);
   dContact contact;
 
   contact.surface.mode = dContactBounce;
   contact.surface.mu = dInfinity;
   contact.surface.bounce = 0.5;
   contact.surface.bounce_vel = 0.001;
 
   if (int num = dCollide(o1,o2,1,&contact.geom,sizeof(dContact))) {
      dJointID c = dJointCreateContact(world,collisions,&contact);
      dJointAttach(c,b1,b2);
   }
 
}

As you can see, instead of doing hard-body collisions between objects, ODE uses 'joints' that simulate a restorative force, pushing the particles back apart. While this is fine for video games and 3D renderings, it falls a bit flat when compared to Dr. Lewis's 'hard-body' collision code. In addition, this joint group has to be created and emptied every timestep, which probably is a nightmare on performance speed. The simloop function handles velocity and position updating in addition to collision handling, and was written to simplify the main simulation loop code.

static void simLoop(int pause) {
   const dReal *pos;
   const dReal *R;
   dSpaceCollide(space,0,&collisionCheck);
   dWorldStep(world,0.000001);
   dJointGroupEmpty(collisions);
   pos = dGeomGetPosition(geom);
   R = dGeomGetRotation(geom);
}

Finally, below is the main loop. It generates 4000 spheres and places them randomly within the bounding box. Note that the spheres it creates are 1e-8 units in size; if you move them to 1e-7 the simulation crashes during the first timestep due to the sheer number of collisions of particles with no velocity. (When the simulation goes to normalize the velocity vector, it fails miserably by dividing by zero and then procuded NaNs for the rest of the runtime.)

// Initialize all of the particles
for(int i = 0; i < 4000; i++) {
      bodies[numBodies] = dBodyCreate(world);
      geom = dCreateSphere(space, 1e-8);  
      dMassSetSphere(&mass,1,1e-7);
      dBodySetMass(bodies[numBodies],&mass);
      dGeomSetBody(geom,bodies[numBodies]);
      dReal x = randPos();
      dReal y = randPos();
      dReal z = randPos()/2;
      dBodySetPosition(bodies[numBodies],x,y,z);
      dBodySetLinearVel(bodies[numBodies],0,-1.5*x,0);
      numBodies++;
   }
 
for(int i = 0; i < numSteps; i++) {
    for(int j = 0; j < numBodies; j++) {
        // Get Position and Velocity
        const dReal *pos = dBodyGetPosition(bodies[j]);
        const dReal *vel = dBodyGetLinearVel (bodies[j]);
 
        // Equations of Motion for Saturn's Gravity
        dReal xAccel = 2*vel[1]+3*pos[0];
        dReal yAccel = -2* vel[0];
        dReal zAccel = -1*pos[2];
 
        // Get Mass for Force computation   
        dBodyGetMass(bodies[j],&mass);
 
        // Compute Forces and Apply
        dBodyAddForce(bodies[j],xAccel*mass.mass,yAccel*mass.mass,zAccel*mass.mass);
 
        if (pos[1] > 1e-6) dBodySetPosition(bodies[j],pos[0],(pos[1]-2e-6),pos[2]);
        if (pos[1] < -1e-6) dBodySetPosition(bodies[j],pos[0],(pos[1]+2e-6),pos[2]);
    }
    // Does Collision and Velocity Updates
    simLoop(i);
}

As is clear, the engine is relatively easy to interact with. Note the simLoop code, which was written to simplify the code for the main loop, and the collision callback function, which must be user-defined in any given code and is called any time there may be a collision between particles.

Polyhedral Simulations

Compilation

Note that getting the ODE code to compile requires multiple library links, and my compilation commands ended up as follows:

g++ -DHAVE_CONFIG_H -I. -I../../ode/src -I/home/cswords/Engine/ode-0.11.1/include -DDRAWSTUFF_TEXTURE_PATH="\"/home/cswords/Engine/ode-0.11.1/drawstuff/textures\"" -DdTRIMESH_ENABLED -DdDOUBLE -g -O2 -MT partSim.o -MD -MP -MF .deps/partSim.Tpo -c -o partSim.o partSim.cpp

mv -f .deps/partSim.Tpo .deps/partSim.Po

/bin/sh ../../libtool —tag=CXX —mode=link g++ -g -O2 -lSM -lICE -L/usr/lib64 -o partSim partSim.o ../../drawstuff/src/libdrawstuff.la ../../ode/src/libode.la -lGLU -lGL -lm -lpthread

g++ -g -O2 -o partSim partSim.o -lSM -lICE -L/usr/lib64 ../../drawstuff/src/.libs/libdrawstuff.a ../../ode/src/.libs/libode.a -lGLU -lGL -lm -lpthread

Problems with ODE

Unfortunately, prolonged simulations run in ODE using basic spheres ran into two major problems: first, the simulations had a tendency to crash after a given amount of time or any time particles were too closely packed. This was a result of vector normalizations dividing by zero and either crashing ODE or producing NaNs for a large number of particle position values. The second, and perhaps more serious, problem with ODE was in its physical accuracy. Wisdom and Tremaine clearly specify expected velocity ellipsoid results which ODE failed to replicate. Given the instability of the engine and the disappointingly inaccurate integrator, we concluded that the Open Dynamcis Engine was not a viable option for polyhedral collision simulations.

Vector Normalization Problems

After any extended run, ODE eventually threw a normalization error and segfaulted to a firy death. Most spherical simulations failed somewhere through the 7th or 8th orbit, and polyhedral simulations failed during the first orbit, thus rendering the polyhedral simulation completely useless; most research simulations performed my Mark Lewis run for between 9 and 12 orbits in order to ensure a stable system before analysing output data.

Velocity Ellipsoid Research using ODE

Velocity ellipsoids for spherical ring simulations are provided in the Wisdom and Tremaine paper, and were used as a heuristic for checking accuracy of the simulation. Unfortunately, as you can see from the image below, the velocity ellipsoids never flatten out, but instead continues to increase or decrease. This implies that ODE's integrator is not scientifically accurate, even with a small enough time step.