Increasing efficiency of N-Body gravity simulation

Question

I'm making a space exploration type game, it will have many planets and other objects that will all have realistic gravity. I currently have a system in place that works, but if the number of planets goes above 70, the FPS decreases an practically exponential rates. I'm making it in C# and XNA.

My guess is that I should be able to do gravity calculations between 100 objects without this kind of strain, so clearly my method is not as efficient as it should be.

I have two files, Gravity.cs and EntityEngine.cs. Gravity manages JUST the gravity calculations, EntityEngine creates an instance of Gravity and runs it, along with other entity related methods.

EntityEngine.cs

        public void Update()
        {
            foreach (KeyValuePair<string, Entity> e in Entities)
            {
                e.Value.Update();
            }

            gravity.Update();
        }

(Only relevant piece of code from EntityEngine, self explanatory. When an instance of Gravity is made in entityEngine, it passes itself (this) into it, so that gravity can have access to entityEngine.Entities (a dictionary of all planet objects))

Gravity.cs

namespace ExplorationEngine
{
    public class Gravity
    {
        private EntityEngine entityEngine;
        private Vector2 Force;
        private Vector2 VecForce;
        private float distance;
        private float mult;

        public Gravity(EntityEngine e)
        {
            entityEngine = e;
        }


        public void Update()
        {
            //First loop
            foreach (KeyValuePair<string, Entity> e in entityEngine.Entities)
            {
            //Reset the force vector
            Force = new Vector2();

                //Second loop
                foreach (KeyValuePair<string, Entity> e2 in entityEngine.Entities)
                {
                    //Make sure the second value is not the current value from the first loop
                    if (e2.Value != e.Value )
                    {
                        //Find the distance between the two objects. Because Fg = G * ((M1 * M2) / r^2), using Vector2.Distance() and then squaring it
                        //is pointless and inefficient because distance uses a sqrt, squaring the result simple cancels that sqrt.
                        distance = Vector2.DistanceSquared(e2.Value.Position, e.Value.Position);

                        //This makes sure that two planets do not attract eachother if they are touching, completely unnecessary when I add collision,
                        //For now it just makes it so that the planets are not glitchy, performance is not significantly improved by removing this IF
                        if (Math.Sqrt(distance) > (e.Value.Texture.Width / 2 + e2.Value.Texture.Width / 2))
                        {
                            //Calculate the magnitude of Fg (I'm using my own gravitational constant (G) for the sake of time (I know it's 1 at the moment, but I've been changing it)
                            mult = 1.0f * ((e.Value.Mass * e2.Value.Mass) / distance);

                            //Calculate the direction of the force, simply subtracting the positions and normalizing works, this fixes diagonal vectors
                            //from having a larger value, and basically makes VecForce a direction.
                            VecForce = e2.Value.Position - e.Value.Position;
                            VecForce.Normalize();

                            //Add the vector for each planet in the second loop to a force var.
                            Force = Vector2.Add(Force, VecForce * mult);
                            //I have tried Force += VecForce * mult, and have not noticed much of an increase in speed.
                        }
                    }
                }

                //Add that force to the first loop's planet's position (later on I'll instead add to acceleration, to account for inertia)
                e.Value.Position += Force;
            }

        }

    }
}

I have used various tips (about gravity optimizing, not threading) from THIS question (that I made yesterday). I've made this gravity method (Gravity.Update) as efficient as I know how to make it. This O(N^2) algorithm still seems to be eating up all of my CPU power though.

Here is a LINK (google drive, go to File > download, keep .Exe with the content folder, you will need XNA Framework 4.0 Redist. if you don't already have it) to the current version of my game. Left click makes a planet, right click removes the last planet. Mouse moves the camera, scroll wheel zooms in and out. Watch the FPS and Planet Count to see what I mean about performance issues past 70 planets. (ALL 70 planets must be moving, I've had 100 stationary planets and only 5 or so moving ones while still having 300 fps, the issue arises when 70+ are moving around)

After 70 planets are made, performance tanks exponentially. With < 70 planets, I get 330 fps (I have it capped at 300). At 90 planets, the FPS is about 2, more than that and it sticks around at 0 FPS. Strangely enough, when all planets are stationary, the FPS climbs back up to around 300, but as soon as something moves, it goes right back down to what it was, I have no systems in place to make this happen, it just does.

I considered multithreading, but that previous question I asked taught me a thing or two, and I see now that that's not a viable option.

I've also thought maybe I could do the calculations on my GPU instead, though I don't think it should be necessary. I also do not know how to do this, it is not a simple concept and I want to avoid it unless someone knows a really noob friendly simple way to do it that will work for an n-body gravity calculation. (I have an NVidia gtx 660)

Lastly I've considered using a quadtree type system. (Barnes Hut simulation) I've been told (in the previous question) that this is a good method that is commonly used, and it seems logical and straightforward, however the implementation is way over my head and I haven't found a good tutorial for C# yet that explains it in a way I can understand, or uses code I can eventually figure out.

So my question is this: How can I make my gravity method more efficient, allowing me to use more than 100 objects (I can render 1000 planets with constant 300+ FPS without gravity calculations), and if I can't do much to improve performance (including some kind of quadtree system), could I use my GPU to do the calculations?

Answer 1

Your approach and implementation are valid (disregarding that position += force line). As long as there is only one instance of gravity. There is nothing in the code you have provided that leads to an over-O(n²) runtime. Therefore it is reasonable to expect a drop in framerate by half at 100 planets. That is 100 planets at 150fps. 200 planets at 30fps. etc.

As this does not match with your observations some questions arise. Most prominently: Are you sure the fps are lost in gravity.Update? Have you measured it with a profiler? Which statement is the offender? Is there any notable difference in memory consumption? Is the garbage collector busy? Can you render 100-200 planets at 300fps when there is no gravity?

Answer 2

You could easily double your speed by adding the gravity force calculated to BOTH objects involved.

You are calculating the force between A & B and the later calculating the force between B & A, but of course, it's the same force in both cases. No need to calculate it twice.

To do this you will need to restructure your loops so that you are calculating the gravity between e1 and only the objects AFTER e1 in your entity list.

//pseudo code
for(i=0 to num_objects){
    for(j=i+1 to num_objects){ //inner loop starts just after position in outer loop 
        f=calc_force_between_objects(i,j)
        i.forces+=f;
        j.forces+=-f; // same force, opposite direction
    }
}

Answer 3

If your universe is divided up in to natural clusters (e.g. a number of solar-systems) Then you can treat each cluster as a single gravitational body for the purposes of gravity simulation (just sum the masses of the objects in the cluster and do a weighted(!) average their positions to get the centre of mass of the cluster). If the cluster is dominated by a large body, e.g. a sun, the just use that as the centre of mass.

This will work fairly accurately for bodies not close the the cluster.

You may also have clusters within clusters. e.g Jupiter and it's moons, is a cluster within the solar system.

Answer 4

Ok I am going to give you simplified answer here. Start with creating a spatial structure like a quadtree. There are plenty of C# docs about this online.

The main idea is to divide your galaxy into a grid, you start with one big cell and divide it into 4 evenly sized cells and so on. Here is an article about QT in C#, but please try to unerstand the concept and don´t copy - paste, as your problem calls for a very special implementation.

After reading the article you should be aware of a mechanism that lets particles register with cells if they move. So if particle A moves from cell [1,1] to cell [2,1], then it has to unregister itself in the first cell and register itself in the second cell.

With this structure you gain a tree, with different granularity. What you can now do is creating a Level of Detail for your mass points. On the finest level (the leafs of the tree) each particle will have its own mass and center. If you now step one level up, you will have a node which contains a number of leafs (planets). You now calculate the center of mass and the mass for all planets combined (by c*alculating the median position of all planets* and adding the mass up).

Now you again step a level up and got a node containing a number of sub nodes (which we just calculated). Now you do the same thing again and generate a new center of mass and a new added mass. Continue with this process until you reach the root node. You will now have a tree, where all planets are singel leafs and every node contains the median center and cummulated mass of all planets below it.

In your Update-logic you now won´t check against every single planet any more. You will calculate the exact force for planets close to the object (maybe the same node level). For planets further away you will not step down into the tree until you reach leafs, instead you will use the previously calculated median values to genereate and approximated force.

Whenever a planet / object changes its tree-cell you need to recalculate the tree from there upwards, but it should still reduce your calculations by a great amound.

Try to implement a system like this and play with the paramters. You should be able to create a quite powerfull system with this tools.

This is a lot of text, try to wrap your head around this and feel free to comment back for questions :)

PS: Its baffeling me that your movings planets create such a huge drop, as you are calculating gravity for static ones too (aren´t you?) maybe there is a also a problem in the movement code.

PS2: maybe try changing if (Math.Sqrt(distance) > (e.Value.Texture.Width / 2 + e2.Value.Texture.Width / 2))

into

if (distance > ((e.Value.Texture.Width * e.Value.Texture.Width ) / 2 + (e2.Value.Texture.Width * e2.Value.Texture.Width) / 2)). while this won´t be the source of your problem it should increase speed a bit ;)

Answer 5

The issue is definitely the N^2 order of complexity. My suggestion is to combine octrees and approximations to avoid as much computation as possible. Some areas where approximation could be used are:

• If m/r^2 is too low, ignore the contribution as too small.
• Any masses within a sufficiently small area, sufficiently far away, can be considered a single mass.

I don't quite have an algorithm yet, but an octree would allow you to ignore a whole node that contributes too little, and also would allow you to treat a concentrated enough node as a single mass without exploring any deeper. "Concentrated enough" could be evaluated in many ways, but a simple formula might just be the area of the node divided by the distance to the node.

Answer 6

This answer is similar to the quad/ oct tree ideas, but not quite the same.

I suspect that there is potential for a massive gain by re-evaluating the force vectors between pairs of bodies based on their proximity.

That is, have each body maintain a priority list of forces applied by each other body. [ This is not really that much storage for a few hundred bodies.]

Having computed the distance, also compute the relative velocity and determine a "good for" time for the force. For example, a force is good as long as its direction does not change by more than 1e-5 radians or its relative distance by more than 1e-5.

class Force
    {
        public Vector2 force;
        public Time    validUntil;
        public Entity  otherBody;
    };

Then in Update, remove those Forces from the head of the priority queue whose validUntil has expired (less than current simulation time). For each of these subtract the force from that currently on the Entity; recompute the force with the otherBody; add that back to the Entity; compute the validUntil and re-insert.

This should vastly reduce the number of times distant objects recompute their force. In fact, you should find that two distant effectively exert constant force on each other in a natural way.

A second issue applies to nearby objects. Use the equivalent of s = ut + at^2/2:

Position oldP;
Velocity oldV;
Acceleration a;
Time dt;

Velocity newV = oldV + dt*a;
Position newApprox = oldP + dt * oldV + 0.5*dt*dt*a;

This will give you much better behavior than your current expression, allowing you to increase dt.

When building a gravity model, I like to apply various tests.

The simplest is to check that you actually get a circular orbit for something like an Earth-Sun system.

More sophisticated is to compute things like total momentum and total angular momentum and verify that these remain sufficiently constant.

Answer 7

The fastest code is code that doesn't run.

Assumption: The realistic gravity is almost surely for objects that are dynamic (spacecraft), where the planets are likely not interacting with each other to cause disturbances (such as our Solar System). Thus, if you have no rogue stars / planets knocking things out of orbit, consider:

Solution: Compute the stable orbit of each moon to its planet, and each planet to its star. You could use Kepler's laws of planetary motion, or even better, assume circular orbits, and just compute the speed required for orbit. Dynamic objects can still use the real gravity formula, so you could recreate the Apollo missions to the Moon, for example.