Like Tree1Likes
  • 1 Post By vart

Gravity simulator implementation questions.

This is a discussion on Gravity simulator implementation questions. within the C++ Programming forums, part of the General Programming Boards category; I've been working on creating a simulator to crash two galaxies together as part of a project to stress test ...

  1. #1
    Registered User TDeiaco's Avatar
    Join Date
    Jul 2013
    Posts
    4

    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.
    Last edited by TDeiaco; 07-09-2013 at 10:23 PM.

  2. #2
    CSharpener vart's Avatar
    Join Date
    Oct 2006
    Location
    Rishon LeZion, Israel
    Posts
    6,484
    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.
    iMalc likes this.
    The first 90% of a project takes 90% of the time,
    the last 10% takes the other 90% of the time.

  3. #3
    Registered User TDeiaco's Avatar
    Join Date
    Jul 2013
    Posts
    4
    Quote Originally Posted by vart View Post
    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;
    
    
        }

  4. #4
    Algorithm Dissector iMalc's Avatar
    Join Date
    Dec 2005
    Location
    New Zealand
    Posts
    6,304
    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.
    My homepage
    Advice: Take only as directed - If symptoms persist, please see your debugger

    Linus Torvalds: "But it clearly is the only right way. The fact that everybody else does it some other way only means that they are wrong"

  5. #5
    - - - - - - - - oogabooga's Avatar
    Join Date
    Jan 2008
    Posts
    2,808
    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.)
    Last edited by oogabooga; 07-10-2013 at 01:22 AM. Reason: not thinking

  6. #6
    Registered User TDeiaco's Avatar
    Join Date
    Jul 2013
    Posts
    4
    Thank you iMalc for your insight and direction I'll give your suggestions a try.

    Quote Originally Posted by oogabooga View Post
    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!!!

  7. #7
    CSharpener vart's Avatar
    Join Date
    Oct 2006
    Location
    Rishon LeZion, Israel
    Posts
    6,484
    I think it should be just sqrt(r) since r is already Distance*Distance
    The first 90% of a project takes 90% of the time,
    the last 10% takes the other 90% of the time.

  8. #8
    Registered User MutantJohn's Avatar
    Join Date
    Feb 2013
    Posts
    1,058
    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...

  9. #9
    - - - - - - - - oogabooga's Avatar
    Join Date
    Jan 2008
    Posts
    2,808
    Quote Originally Posted by TDeiaco View Post
    Can you explain why you sqrt(r*r); because wouldn't this simple produce 'r' again?
    Yes. It should be sqrt(r).

  10. #10
    Registered User TDeiaco's Avatar
    Join Date
    Jul 2013
    Posts
    4
    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.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. SDL gravity.
    By bijan311 in forum Tech Board
    Replies: 4
    Last Post: 06-02-2010, 08:08 PM
  2. questions about classes(implementation files)
    By Link_26 in forum C++ Programming
    Replies: 6
    Last Post: 06-22-2006, 09:48 AM
  3. Need help with gravity engine
    By madmech in forum Game Programming
    Replies: 11
    Last Post: 07-03-2005, 02:41 PM
  4. Gravity pong!
    By Nutshell in forum Game Programming
    Replies: 10
    Last Post: 04-25-2003, 08:52 AM
  5. Gravity? We don't need no steenking gravity!
    By Govtcheez in forum A Brief History of Cprogramming.com
    Replies: 22
    Last Post: 03-28-2002, 06:24 PM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21