Thread: Gravity simulator implementation questions.

I've been working on creating a simulator to crash two galaxies together as part of a project to stress test a CUDA super computer. I've got a long way to go and am currently just working on correctly simulating n-body gravity functions. First I will use this to simulate the cores of the galaxies (the black holes) and eventually the stars.

So long story short I'm working on the beginnings of a gravity simulator. At this point I found some basic code that works well but doesn't quite give the effect I'm looking for.

The code below only pulls each object towards each other like a spring faster and faster until they shoot off into infinity. I try to give one of my bodies an initial velocity to get it to orbit another, but it always just shoots straight at the other body. I'm thinking I need to factor in inertia so that the initial velocity doesn't just get calculated away really fast by the other calculations.

I'm really looking for a bit of direction to get a real gravity simulator with orbits and such working right so eventually I can scale it up to a galaxy, throw in 100B stars and let the CUDA run for a month..

Code:

void update_galaxies(GLdouble elapsedTime)    {
        //Calculate gravity simulations
        GLdouble r1, r2, r3;
        r1 = r2 = r3 = 0.0;


        for(unsigned int i = 0; i < galaxies.size(); i++)
        {
            for(unsigned int j = 0; j < galaxies.size(); j++)
            {
                //Don't caculate forces on the same galaxies
                if(i != j){
                    r1 = galaxies[i]->xpos - galaxies[j]->xpos; 
                    r2 = galaxies[i]->ypos - galaxies[j]->ypos;
                    r3 = galaxies[i]->zpos - galaxies[j]->zpos;


                    galaxies[i]->xacc += -r1 * GRAV * galaxies[j]->mass / ((r1*r1 + r2*r2 + r3*r3)/3.0)*((r1*r1 + r2*r2 + r3*r3)/3.0);
                    galaxies[i]->yacc += -r2 * GRAV * galaxies[j]->mass / ((r1*r1 + r2*r2 + r3*r3)/3.0)*((r1*r1 + r2*r2 + r3*r3)/3.0);
                    galaxies[i]->zacc += -r3 * GRAV * galaxies[j]->mass / ((r1*r1 + r2*r2 + r3*r3)/3.0)*((r1*r1 + r2*r2 + r3*r3)/3.0);
                }
            }
        }


        for(int i = 0; i < galaxies.size(); i++)
        {
            galaxies[i]->update_pos(elapsedTime);
        }
    }

As you can see, I'm calculating all the bodies in a vector called "galaxies" with each other, and doing a basic gravity calculation to it. The update_position function simply takes the calculated acceleration and uses it to calculate the velocity and position based on the "elapsedTime". Pretty simple but ask any questions that would clear up what I'm doing.

I think I need to use the Varlet or Runge-Kutta integration methods, after doing a bit more research. A bit of direction would still be great, but I think after some research on these integration methods I will have something a bit more accurate.

I do not see in your code velocity vector at all

Gravity used to calculate acceleration.
Acceleration used to calculate change in velocity vector
Velocity vector is used to calculate change in coordinates.

That's how I would do it.

No it's not inertia you need. If the initial velocity of one of them is tangential to the vector between the objects, then so long as you haven't added the effect of friction into the simulation, then one object will simply orbit the other. If that initial velocity is too low or too high then that orbit will just be eliptical.

Shooting straight into the other body could just mean that initial velocity is so low that the eliptical orbit becomes almost flattened to a line, or the mass is too large. Try a larger initial velocity, or a smaller mass.
Make sure that things are in compatible units E.g. don't mix meter's per second per second and kilometres per hour.

Also lines 18-20 are certainly incorrect as the later bracketed parts cancel out and thus they are equivalent to this:

Code:

                    galaxies[i]->xacc += -r1 * GRAV * galaxies[j]->mass;
                    galaxies[i]->yacc += -r2 * GRAV * galaxies[j]->mass;
                    galaxies[i]->zacc += -r3 * GRAV * galaxies[j]->mass;

Consider using a mathematical vector made of three GLdouble's that have operator overloads for addition, and methods for length etc, so that you can do calculations such as the above mostly on one line. Or at the very least, factor out the common length calculation such that it is only done once per inner loop iteration, rather than relying on the compiler to do that for you. It greatly increases readability anyway.

This seems wrong:

Code:

galaxies[i]->xacc += -r1 * GRAV * galaxies[j]->mass / ((r1*r1 + r2*r2 + r3*r3)/3.0)*((r1*r1 + r2*r2 + r3*r3)/3.0);

Whats's up with dividing by 3? Where is the i galaxy's mass? Maybe something clever is going on there that I don't understand, but why don't you try a more straightforward method like this:

Code:

r = (r1*r1+r2*r2+r3*r3);
sr = sqrt(r*r);
a = GRAV * galaxies[i]->mass * galaxies[j]->mass / r;
galaxies[i]->xacc += a * r1 / sr;
galaxies[i]->yacc += a * r2 / sr;
galaxies[i]->zacc += a * r3 / sr;

(It's also possible that you need to reverse the signs on r1,r2,r3.)

I think it should be just sqrt(r) since r is already Distance*Distance