# Optimizing gravity calculations

I've got a bunch of objects of varying size and velocity which gravitate towards each other. On every update, I have to go over every object and add up the forces due to gravity of every other object. It doesn't scale very well, is one of two big bottlenecks I've found in my game, and I'm not sure what to do to improve the performance.

It feels like I should be able to improve the performance. At any given time, probably 99% of the objects in the system will only have a negligible influence on an object. I of course can't sort the objects by mass and only consider the top 10 largest objects or something, because the force varies with the distance more than with the mass (the equation is along the lines of force = mass1 * mass2 / distance^2). I think a good approximation would be to consider the largest objects and the closest objects, ignoring the hundreds of tiny fragments of rock on the other side of the world that can't possibly affect anything -- but in order to find out which objects are closest I have to iterate over all the objects, and their positions are changing constantly, so it's not like I can just do it once.

Currently I'm doing something like this:

private void UpdateBodies(List<GravitatingObject> bodies, GameTime gameTime)
{
for (int i = 0; i < bodies.Count; i++)
{
bodies[i].Update(i);
}
}

//...

public virtual void Update(int systemIndex)
{
for (int i = systemIndex + 1; i < system.MassiveBodies.Count; i++)
{
GravitatingObject body = system.MassiveBodies[i];

Vector2 force = Gravity.ForceUnderGravity(body, this);
ForceOfGravity += force;
body.ForceOfGravity += -force;
}

Vector2 acceleration = Motion.Acceleration(ForceOfGravity, Mass);
ForceOfGravity = Vector2.Zero;

Velocity += Motion.Velocity(acceleration, elapsedTime);
Position += Motion.Position(Velocity, elapsedTime);
}


(note that I've removed a lot of code -- for example the collision tests, I do not iterate over the objects a second time to detect collisions).

So I'm not always iterating over the whole list -- I only do that for the first object, and every time the object finds the force it feels toward another object, that other object feels the same force, so it just updates both of them -- and then that first object doesn't have to be considered again for the rest of the update.

The Gravity.ForceUnderGravity(...) and Motion.Velocity(...), etc. functions just use a bit of XNA's built in vector math.

When two objects collide, they create massless debris. It's kept in a separate list and the massive objects do not iterate over the debris as part of their velocity calculation, but each piece of debris must iterate over the massive particles.

This doesn't have to scale to incredible limits. The world isn't unlimited, it contains a border which destroys objects that cross it -- I'd like to be able to handle maybe a thousand or so objects, currently the game starts to choke around 200.

Any thoughts on how I might improve this? Some heuristic I can use to shave the loop length from hundreds down to just a few? Some code I can execute less often than every update? Should I just multithread it until it's fast enough to allow for a decent sized world? Should I try to offload the velocity calculations to the GPU? If so, how would I architect that? Can I keep static, shared data on the GPU? Can I create HLSL functions on the GPU and call them arbitrarily (using XNA) or do they have to be part of the draw process?

-
Just a note, you said "the massive objects do not iterate over the debris as part of their velocity calculation, but each piece of debris must iterate over the massive particles." I got the impression from that, that you're assuming it's more efficient. However, iterating 100 debris objects 10 times each, is still the same as iterating 10 massive objects 100 times. Perhaps, it would be a good idea to iterate each debris object in the massive objects loop so you don't do it a second time. –  Richard Marskell - Drackir Nov 4 '11 at 19:05
How accurate a simulation do you need? Do you really need everything accelerating toward each other? And do you really need to use a true gravitation calculation? Or can you depart from those conventions for what you're trying to accomplish? –  chaosTechnician Nov 4 '11 at 19:09
@Drackir I think you're right. Part of the reason they're separated is because the math is different for massless objects, and part is because they originally didn't obey gravity at all, so it was more efficient to not include them. So it's vestigial –  Carson Myers Nov 4 '11 at 19:15
@chaosTechnician it doesn't have to be very accurate -- in fact if it's only accounting for the few most dominant forces then a system would be more stable, which is ideal. But it's finding out which of the forces are most dominant in an efficient way that I'm having trouble with. Also the gravitation calculation is already approximated, it's just G * m1 * m2 / r^2, where G is just to tweak the behavior. (although I can't just have them following a path, because the user can disturb the system) –  Carson Myers Nov 4 '11 at 19:16
Why does each piece of debris have to iterate over the massive particles if it's massless? Collisions? –  Sam Hocevar Nov 4 '11 at 20:18

This sounds like a job for a grid. Divide your game space into a grid and for each grid cell keep a list of the objects currently in it. When objects move across a cell boundary, update which list they're in. When updating an object and searching for others to interact with, you can look at just the current grid cell and a few neighboring ones. You can tweak the size of the grid for the best performance (balancing the cost of updating the grid cells - which is higher when the grid cells are too small - with the cost of doing the searches, which is higher when the grid cells are too large).

This will, of course, cause objects that are farther apart than a couple of grid cells not to interact at all, which is probably an issue because a large accumulation of mass (either a big object, or a cluster of many small objects) should, as you mentioned, have a larger region of influence.

One thing you could do is keep track of the total mass within each grid cell, and treat the whole cell as a single object for the purposes of farther-away interactions. That is: when you calculate the force on an object, calculate the direct object-to-object acceleration for the objects in a few nearby grid cells, then add in a cell-to-cell acceleration for each farther-away grid cell (or maybe just the ones with a non-negligible amount of mass in them). By a cell-to-cell acceleration, I mean a vector computed using the two cells' total masses and the distance between their centers. That should give a reasonable approximation of the summed gravity from all the objects in that grid cell, but much more cheaply.

If the game world is very large, you could even use a hierarchical grid, like a quadtree (2D) or octree (3D), and apply similiar principles. Longer-distance interactions would correspond to higher levels of the hierarchy.

-
+1 for the grid idea. I'd suggest also tracking the center of mass for the grid as well to keep the calculations a little more pure (if necessary). –  chaosTechnician Nov 4 '11 at 19:13
I like this idea a lot. I had considered containing the objects in cells, but abandoned it when considering two nearby objects which were technically in different cells -- but I didn't make the mental leap to consider a few adjacent cells, as well as considering the combined mass of other cells. I think this should work really well if I do it right. –  Carson Myers Nov 4 '11 at 19:30
This is essentially the Barnes-Hut algorithm: en.wikipedia.org/wiki/Barnes–Hut_simulation –  Russell Borogove Nov 4 '11 at 20:24
It even sounds similar to how gravity is supposed to work in physics -- the bending of space-time. –  Zan Lynx Nov 4 '11 at 20:52
I like the idea - but I honestly think it needs a tiny bit more refinement - what happens if two objects are very close to each other but in separate cells? I could see how your third paragraph might help - but why not just do a circular cull on the interacting objects? –  Jonathan Dickinson Jan 12 '12 at 17:01

The Barnes-Hut algorithm is the way to go with this one. It's been used in supercomputer simulations to solve your exact problem. It's not too hard to code, and it's very efficient. I actually wrote a Java applet not too long ago to solve this problem.

Visit http://mathandcode.com/programs/javagrav/ and press "start" and "show quadtree".

In the options tab, you can see that the particle count can go all the way up to 200,000. On my computer, the calculation finishes in about 2 seconds (the drawing of 200,000 dots takes about 1 second, but the calculation runs on a separate thread).

Here's how my applet works:

1. Create a list of random particles, with random masses, positions, and starting velocities.
2. Construct a quadtree from these particles. Each quadtree node contains the center of mass of each subnode. Basically, for each node you have three values: massx, massy, and mass. Every time you add a particle to a given node, you increase massx and massy by particle.x*particle.mass and particle.y*particle.mass respectively. The position (massx/mass, massy/mass) will end up as the node's center of mass.
3. For each particle, calculate the forces (fully described here). This is done by starting at the top node and recursing through each subnode of the quadtree until the given subnode is small enough. Once you stop recursing, you can calculate the distance from the particle to the node's center of mass, and then you can calculate the force using the node's mass and the particle's mass.

Your game should easily be able to handle a thousand mutually attracting objects. If each object is "dumb" (like the bare-bones particles in my applet), you should be able to get 8000 to 9000 particles, maybe more. And this is assuming single-threading. With multi-threaded or parallel computing applications, you can get many more particles than that updating in real time.

-
The first link is dead. Is the applet hosted elsewhere? –  Anko Feb 9 '13 at 11:15
Fixed! sorry, forgot to pay the rent on that domain, and someone auto-bought it :\ Also, 3 mins is a pretty good response time on a 1.3 year old post 8D –  NeuroFuzzy Feb 9 '13 at 11:20
And I should add: I don't have the source code. If you're looking for some source code check out part-nd (written in c). I'm sure there are others out there, too. –  NeuroFuzzy Feb 9 '13 at 11:22

Nathan Reed has an excellent answer. The short version of it is to use a broadphase technique that fits your simulation's topology, and to only run the gravity calculations on pairs of objects that are going to have a noticeable effect on each other. It's really no different than what you'd do for regular collision detection broadphase.

Carrying on from that, though, another possibility is to only update objects intermittently. Basically, each time step (frame) only update a fraction of all the objects, and leave the velocity (or acceleration, depending on your preference) the same for the other objects. The user is unlikely to notice any delay in the updates from this so long as the intervals aren't too long. This will give you a linear speed up of the algorithm, so definitely look into the broadphase techniques like Nathan suggested as well, which can give much more significant speedups if you have a ton of objects. While not at all modeled the same in the slightest, this is kind of like having "gravity waves." :)

Also also, you could generate a gravity field in one pass, then update the objects in a second pass. The first pass you're basically filling a grid (or a more complex spatial data structure) with the gravity influences of each object. The result now is a gravity field that you can even render (looks pretty cool) to see what acceleration will be applied to an object at any given location. Then you iterate over the objects and simply apply the gravity field's effects to that object. Even cooler, you can do this on a GPU by rendering the objects as circles/spheres to a texture, then reading the texture (or using another transform-feedback pass on the GPU) to modify the object's velocities. So long as the GPU has a decent fillrate and your objects have a reasonable (if technically unrealistic) cut-off distance for their gravitational effect, this will likely be the fastest approach (but profile and test, as always).

-
Separating the process into passes is an excellent idea, since the update interval is (as far as I know) a very small fraction of a second. The gravity field texture is AWESOME but perhaps a little beyond my reach right now. –  Carson Myers Nov 4 '11 at 19:37
Don't forget to multiply the applied forces by how many time-slices have been skipped since the last update. –  Zan Lynx Nov 4 '11 at 20:53
@seanmiddlemitch: Could you elaborate on the gravity field texture a little more? I'm sorry, I'm not a graphics programmer, but this sounds really interesting; I just don't understand how it's supposed to work. And/or maybe you have a link to a description of the technique? –  Felix Dombek Nov 5 '11 at 21:10
@FelixDombek: Render your objects as circles representing area of influence. Fragment shader writes a vector pointing at the center of the object and with appropriate magnitude (based on distance from center and mass of object). Hardware blending can handle summing these vectors in additive mode. The result will not be accurate gravity fields, but will almost certainly be good enough for a game's needs. As yet another approach using the GPU, see this CUDA-based N-body gravity simulation technique: http.developer.nvidia.com/GPUGems3/gpugems3_ch31.html –  Sean Middleditch Nov 15 '11 at 0:10

I would recommend using a Quad Tree. They allow you quickly and efficiently look up all objects in an arbitrary rectangular area. Here's the wiki article on them: http://en.wikipedia.org/wiki/Quadtree

I would also maintain a list of all the large objects so that they can interact with everything regardless of distance.

-
Thanks! I may use this :) –  Carson Myers Nov 4 '11 at 21:50

Just a small bit of (possibly naive) input. I don't do game programming, but what I am feeling is that your fundamental bottleneck is the gravity-due-to-gravity calculation. Instead of iterating over each object X and then finding the gravitational effect from each object Y and adding it, you can take each pair X,Y and find the force between them. That should cut the number of gravity calculations from O(n^2). Then you will be doing a lot of addition (O(n^2)), but this is normally less expensive.

Also at this point you can implement rules such as "if the gravitational force will be less than \epsilon because these bodies are too small, set the force to zero". It may be advantageous to have this structure for other purposes too (including collision detection).

-
This is fundamentally what I'm doing. After I get all pairs involving X, I don't iterate over X again. I find the force between X and Y, X and Z, etc, and apply that force to both objects in the pair. After the loop completes, the ForceOfGravity vector is the sum of all the forces, and that is then converted into a velocity and new position. I'm not sure that the gravity calculation is particularly expensive, and checking if it exceeds a threshold first wouldn't save a noticeable amount of time, I don't think –  Carson Myers Nov 4 '11 at 22:44

In extending seanmiddleditch's answer, I thought I might shed some light (irony?) on the gravity field idea.

Firstly, don't think of it as a texture, but a discrete field of values that can modified (a two-dimensional array, as it were); and the subsequent accuracy of the simulation could be the resolution of that field.

When you introduce an object into the field, its gravitation potential can be calculated for all surrounding values; thereby creating a gravitation sink in the field.

But how many of these points should you calculate before it becomes more or as ineffective as before? Probably not many, even 32x32 is a substantial field to iterate for each object. Therefore break the entire process into multiple passes; each with varying resolutions (or accuracy).

Ie, the first pass may calculate the objects gravity represented in a 4x4 grid, with each cell value representing a 2D coordinate in space. Giving an O(n * 4 * 4) sub-total complexity.

The second pass may more accurate, with a 64x64 resolution gravity field, with each cell value representing a 2D coordinate in space. However, as the complexity is very high, you can restrict the radius of surrounding cells affected (perhaps, only the surrounding 5x5 cells are updated).

An additional third pass could be used for high accuracy calculations, with maybe a resolution of 1024x1024. Remembering at no time are you actually performing 1024x1024 separate calculations, but operating only on portions of this field (perhaps 6x6 sub-sections).

In that way, your overall complexity for the update is O(n * (4*4 + 5*5 + 6*6)).

To then calculate the velocity changes to each of your objects, for each gravity field (4x4, 64x64, 1024x1024) you just map the point masses position to a grid cell, apply that grid cells overall gravitational potential vector to a new vector; repeat for each "layer" or "pass"; then add them together. This should give you a good resultant gravitational force vector.

Therefore, the overall complexity is: O(n * (4*4 + 5*5 + 6*6) + n). What really counts (for complexity) is how many surrounding cells you update when calculating the gravitational potential in the passes, not the overall resolution of the gravity fields.

The reason for low resolution fields (first passes) is to obviously encompass the universe as a whole, and ensure outlying masses are attracted to more dense areas despite distance. Then use higher resolution fields as separate layers to increase accuracy for neighboring planets.

-