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.

Quote:

---

Originally Posted by vart

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.

---

I said I do this in the "update_pos()" function but I'll show it here anyways.

Code:

---

    void update_pos(GLdouble deltaTime){       
        xvel += xacc * deltaTime;
        yvel += yacc * deltaTime;
        zvel += zacc * deltaTime;


        xpos += xvel * deltaTime;
        ypos += yvel * deltaTime;
        zpos += zvel * deltaTime;


    }

---

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.)

Thank you iMalc for your insight and direction I'll give your suggestions a try.

Quote:

---

Originally Posted by oogabooga

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.)

---

Thanks oogaabooga, I was
copying an algorithm from a website to just get it working. Your approach seems much more straightforward and understandable. Can you explain why you sqrt(r*r); because wouldn't this simple produce 'r' again? I know there's a reason and I would like to understand it better.

BTW I used your guys' advice and I'm getting great looking arcs and orbits! Thanks!!!

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

I know this is intended as a stress test but super computers deserve better algorithms than the (n/2)(n-1) calculations required.

If you're interested because I was at the time, you should look up how simulation codes such as Gadget or Enzo work. This is how much larger-scale simulations are performed and will make you much better at computer science if you weren't already good at it.

What they do is, they create a 3D distribution of points in space. They allow multiple particles to reside in the same point but they really split the particle itself into multiple parts spread over the points in space. This is because more often than not a particle is in-between two points in space so it's a better approximation to re-distribute the mass based off of that.

But what's the point of putting particles on a 3D mesh? Simple. The Fastest Fourier Transform in the West. You can solve for gravity with just a density function and FFTW and it's so much faster. In fact, the calculation doesn't even have a derivative, it's literally just multiplying by a constant created by the Fourier transform.

So codes like Enzo will use this but then choose to refine the mesh with separately allocated side-meshes that are solved independenly while something like Gadget will instead create an octree and solve for gravity on that instead of refining a mesh. Such codes are called hybrid codes as the original octree implementation was purely just an octree, originally called a Barnes-Hut tree.

This was all written assuming you're coming from a physics background like me. If not then learn more physics! It makes your calculations better. The graphics programmers benefit from it and shoot, we benefit from the graphics programmers too. Now if only I could get hired...

Quote:

---

Originally Posted by TDeiaco

Can you explain why you sqrt(r*r); because wouldn't this simple produce 'r' again?

---

Yes. It should be sqrt(r).

Thanks Mutant John, that's just the direction I was looking for. I went to school for Electrical Engineering so I really don't have a great background in CS and an even more hobbyist background in Physics. Your advice is great! Thanks for your time.