What's so CPU-intensive about n-body physics?

[JavaProphet]

I've read some posts here about adding n-body physics to the game, and people seem to believe that it's not algorithmically complex, but very CPU intensive? Where does this CPU intensiveness come from? As a hobby programmer, assuming you have your calculation for gravity for any given celestial body, for any given craft(in this case), would you not just calculate the gravity for the other bodies, and add the vectors, to get a final gravity vector? I understand you'd be getting into some high decimal point calculations(easily thousands of zeroes from my first guess), but there are a lot of optimization techniques to be had.

First thing that comes to mind:

Cache the values, radially, for a certain distance from each body, exponentially increasing to infinity(for efficiency with inverse-square). For instance 0-10000 has 9.81g(kerbin), 10000-100000 has 9.8g(example), so on and so forth. A one-time calculation for each distance. Since it's a game, we don't need insane to the power of -10000 accuracy.

I might just be a moron, but what is computationally slow about any of this?

For those interested, the original thread another user posted in the Suggestions section was closed, so I didn't reply to it.

EDIT: I did more research, and people also say on-rails is impossible for this. From what I know, KSP does on-rails by getting position based on time, rather than frame-by-frame. If you could do subframes, or have your gameframe, then a simulation-frame, in which you did a pure lookup for the gravity for each body on your object, added, and then calculated a force vector, treated everything(crafts) as points rather than complex crafts(part-part physics), I don't see so much CPU lag involved. This might require a rework of savegames to instead specify current location and velocity rather than orbit and time, although I personally wouldn't expect/want this in the game for gameplay reasons.

Edited August 20, 2014 by JavaProphet

[Mr Shifty]

1) KSP uses patched conics because it's a good 1st order approximation, conceptually easy to understand, and can be solved analytically (thus projected for arbitrary time-step intervals.) It's not because it's computationally difficult. Orbiter has been doing n-body physics in its spaceflight simulation for over a decade, and HarvesteR played Orbiter long before KSP development was done. He's known what's possible and made a deliberate decision to do it the way it's done.

2) You can't cache gravity values (or more accurately, it doesn't save you any work.) N-body gravitation in a fixed universe (where the planets and stars just hung at fixed positions in space) would be simple to calculate. Trajectories would be fixed, and you could indeed pre-calculate and cache a grid of gravity vector forces. (See former KSP dev NovaSilisko's WIP game [thread=67691]Science of the Spheres[/thread] for an example of this kind of universe.) But in real planetary systems, everything is in constant motion. The gravity gradient is in constant flux, and trajectories at or near Hill sphere boundaries can become pretty complex.

3) The difficult part of n-body physics (especially for KSP) is not calculating where my vessel will go in the next frame; it's projecting it out over a large segment of time in order to draw trajectories on the map view and allow time acceleration over arbitrary time-steps. There's a detailed discussion, and a project working on this, [thread=68502]here[/thread].

[JavaProphet]

1) I completely agree with this. I'm not saying I want n-body physics, but rather curious about it.

2) I'm talk about a per-body cache. May not be necessary, depends on the complexity of the gravity calculation. I'll assume I'm missing some important fact about gravity.

3) For on-rails simulation, you could just have the sim-frame thread go ahead of time to calculate this out. For non-on-rails, you would use the same sim-thread. You could reduce the poll-density as the trajectory grew longer.

I'll take a look at that thread, thanks.

[Kermunist]

Well, first, everything is relative. So Lets have a think about patched conics, what are they like to calculate?

Very easy. Patched conics are analytically solvable for the position with time. That means that given the starting position and velocity of the spacecraft at time 0, I can write the position and velocity at some future time t using a formula. For now I will use R to represent the canonical co-ordinates of the problem, it is a vector like (x,y,z,v_x,v_y,v_z), so I can write:

R = F(R_t=0,t)

Where F can be found in advance and put into the code. This is great, because you can run the simulation at any speed with pretty good accuracy. If you want each time step to be 0.1 seconds (as I believe is the case at timewarp 1) then you take the R from the last iteration, use t=0.1 an do the calculation. It costs some number of operations, lets say N flops. If on the other hand the timewarp is set to 100,000, you do the same calculation but with t=10,000. It still costs N flops, and in both cases the solution is exact to up to machine precision.

Strictly speaking this assumes no SOI changes. Change SOI and you need to change F to the one governing the new SOI. If you just change F for the first step after entering the new SOI, then there is no additional computational cost. I think that is what KSP does, and why it is sometimes inadvisable to cross SOI boundaries at speed. More accurate options, with a higher computational cost, are possible.

Now, there is no such solution for an n-body equation. This is not a computational problem, it is a mathematical one. The mathematicians simply have not found a way of obtaining F for more than two bodies, and a few special cases with 3 bodies. So now you have to do the calculation in numerical steps. So the question is how many steps per in-simulation second are needed to be accurate? If you don't do enough, then small fast orbits get turned into short joined up straight lines, which don't really represent the orbit perfectly. Since each step epends on the one before, these errors add up to a real mess. My gut feeling is that 10 steps per in-game second will be required, and each one will require some number of operations N, probably not so different from the N above. So no change in computational load at 1x warp. But now, when you timewarp, you can't just change t. Every step depends on the output of the step before, and all must be done. So you have to do 100,000 times as many calcuations. This is where it starts to get messy.

For some rough figures pulled out of thin air: Lets assume that there are 50 ships to track, each of which needs calculations for the nearest 3 bodies. All at max timewarp. That requires 50x3x100000=15M calculations per second, instead of 50 needed under patched conics.

Is this impossible? well, maybe. It depends on the calculation being performed. My gut feeling is that it will be beyond most processors. If it is possible it will probably only really work well with some very well optimized code, making careful use of the processor's features. The lookup table you describe for example is probably not good, as it will be bigger than the processor's cache, and have to live in RAM. Fetching data from RAM takes time, direct calculations will probably be faster in this case. I think Orbiter does try to do n-body calculations in this way, but I think it has problems at high warp (I've never played it)

Also, I don't think it would add to the game, as n-body systems don't usually have stable orbits. But this is a question in the science labs about computational complexity, so I won't go into that.

[Mr Shifty]

2) I'm talk about a per-body cache. May not be necessary, depends on the complexity of the gravity calculation. I'll assume I'm missing some important fact about gravity.

I don't think you're missing anything; I suspect that it's just far more difficult to do vector addition (required for summing gravity components for all bodies) than it is to do scalar multiplication (which is all that's required to calculate gravity for a single body.) So the vector stuff dominates, and can't be cached since a vessel's relative position changes constantly. Also, and this is just speculation since I know next to nothing about numeric methods, it's possible that the integrations necessary to do n-body projections require the analytic single-body solutions.

3) For on-rails simulation, you could just have the sim-frame thread go ahead of time to calculate this out. For non-on-rails, you would use the same sim-thread. You could reduce the poll-density as the trajectory grew longer.

I gather that the central issue is how much precision you retain when you reduce the poll density. Conic trajectories can be stored using six scalar parameters. How would you store a Newtonian trajectory and retain enough precision at high time-steps when near a planet to avoid crashing into the surface? (Note this happens at high time-steps in Orbiter.) Most of the exciting part of Principia for me is how it could generalize trajectories, allowing for trajectory generation for ships under constant thrust, aerobraking trajectories, trajectories in rotating (non-galactic) frames, non-spherical gravity, etc. (Note -- and I'm sorry to belabor this point -- all of this is already possible in Orbiter.)

Edited August 20, 2014 by Mr Shifty

[Kermunist]

I think actually the vector addition is not too bad. It is after all just three addition operations.

Depending on the processor design and my ability to dredge these things out of the back of my memory:

- Addition costs about 1 cycle

- Multiplication costs about 3 cycles

- Division costs much more and can vary a lot, 10-30 cycles

So it's probably going to be the division by distance squared which dominates the calculation. Still faster than hitting RAM though, I think.

[Mr Shifty]

I think actually the vector addition is not too bad. It is after all just three addition operations.

Depending on the processor design and my ability to dredge these things out of the back of my memory:

- Addition costs about 1 cycle

- Multiplication costs about 3 cycles

- Division costs much more and can vary a lot, 10-30 cycles

So it's probably going to be the division by distance squared which dominates the calculation. Still faster than hitting RAM though, I think.

How are you doing vector ops with just multiply and add? You're going to need either trigonometric functions or division. Consider r_s is the vector position of my vessel and r_p is the vector position of a planet. You're given the magnitude of the gravitation force between them (or can calculate it easily.) How would you calculate the force vector on my vessel without trig or division?

[JavaProphet]

I made a 2D physics engine, and with it my own vector math library. Vector addition takes the number of computations for scalar addition * number of dimensions. You would only reduce poll density for the trajectory, otherwise you need the same poll density, or your going to have saving/loading issues(if you keep up the abstract time step idea). If you didn't base where the craft is based on time, but rather on a static coordinate plane(which is a problem in itself), and a velocity vector, you wouldn't have to worry about save/load sync. Some calculation math for 100000x warp:

Assume we need 50(ships)*n(in this example, 10)*100000(warp) calculations.

As for a 'compression', you could reduce poll density for time warp, but this would either need to be recorded, or you would need to discard the time step thing.

Each calculation would be calculating gravity (mass*distance*grav_contstant, or 6 cycles according to Kermunist) = 6 cycles.

Accessing an array for an arbitrary entry, would take an average of 5 cycles per access. 5 * 10 = 50 cycles (for cache)

So, my bet is that real-time-calc is better for anything with more than 3 bodies. (1.5 * 3 = 4.5 cycles)

For n = number of planets:

So, 6 cycles per body to get the gravity speed. 6n cycles.

Then, you need to subtract any vessels position from all planets positions, and then multiply by that gravity for each body(to get the gravity vector).

So, 3n cycles + 9n(? multiplication of vectors, not going to think about too critically right now) cycles = 12n cycles to get the gravity vector.

So now we have 18n cycles.

NOTE* at this point, Kermunist's post about the CPU cycles for these ops may not be too accurate for precision floating point ops. We may use a cache or perhaps just more cycles than I'm calculating. I don't work with processor-level stuff much.

Now we need to add up all of those vectors, so 3n more cycles to get the final vector.

A total of 21 * number of planet CPU cycles per ship per calculation.

So for 50 ships, 10 planets, and 100000x warp, we would need (assuming 1 step/calculation/warp) 50*10*100000*21*10 CPU cycles per second, or: 10.5 billion CPU cycles, 10.5 GHz of processor power.

Before you say, we don't have much past 4, we all run an average of 6 cores(or so), so we have 24 GHz available(potentially).

It's not impossibly on your typical gaming computer, but it'd be considerable. However, that is 100000x warp. This does not take into account how we would be dealing with floating point ops, so it could be much higher.

EDIT: I do now see you would need some trig in order to get the actual force vector + you need to apply it as a force and calculate velocity. The real number is probably 4 times higher or more.

Edited August 20, 2014 by JavaProphet

[RockoDyne]

At 100,000x time warp, patched conics only has one set of calculations to perform on all vessels on rails.

Now for n-body physics to be accurate, how many calculations need to be performed per vessel per frame at 100,000x warp? I'm going to guess and say it will be enough to qualify as a computationally expensive operation.

[Kermunist]

How are you doing vector ops with just multiply and add? You're going to need either trigonometric functions or division. Consider r_s is the vector position of my vessel and r_p is the vector position of a planet. You're given the magnitude of the gravitation force between them (or can calculate it easily.) How would you calculate the force vector on my vessel without trig or division?

I was assuming that one would have the vector distance to the planet, use this to calculate the vector force, and then you just add the vector forces. In which case I think you only need three divisions per planet in total (though I haven't sketched it out completely, it might be four), all in calculating the vector force. If aiming for speed, you'd want to stay away from trig, as it's rather costly to evaluate.

edit: I think I can get it with one square root, three divisions, ten multiplications and five additions per planet (to update V, another 3 multiplications and 3 additions or so for R)

@JavaProphet

I'm not sure I follow your post, I may have a closer look later. In reply to your note, the numbers I gave are the typical times it takes for the ALU to complete the requested operation. At the very least, there will be one extra cycle to read the instruction, fetch the operands and pass them to the ALU, and that assumes they are available in the processor. To avoid having any other overheads, and spend every cycle doing what you want, you have to program very carefully. Assembly would be best, compiled languages like C or FORTRAN may also let you get there. I would be amazed if you were anywhere near that efficient in Java (guessing from your name) or C# (which KSP uses).

edit: here's a comparison of the time taken on some operations by different processors http://www.agner.org/optimize/instruction_tables.pdf

Edited August 20, 2014 by Kermunist

[TythosEternal]

Numerical stability is also an issue. With restricted two-body, there's really only one numerical problem to solve at each time-step (anomaly), and there are very efficient algorithms for doing so with a controllable amount of error. When solving an n-body by straight integration (even a restricted two-body, actually, without the right conservative integrators), small errors (round-off in particular) will be exponentially magnified over time.

Computing the EoMs are actually tricky in and of themselves, even ignoring the integration step, because of a) the way your problem scales (though this can be partly addressed by bounding knowledge of your problem size and by using dipoles for influence grouping), that nasty square root operator with which you are evaluating the magnitude of your relative range vectors, and c) the loss of precision you encounter when you sum large numbers of near-uniformly-distributed vectors, many of which will have at least one component opposed to it's brethren.

[Mr Shifty]

I was assuming that one would have the vector distance to the planet, use this to calculate the vector force, and then you just add the vector forces.

Calculating the relative position vectors can be done by subtracting absolute position vectors. That's great, but the only thing the force vector shares with the position vector is direction. In order to apply the force in the correct direction, you either need to use trigonometry or you can divide the relative position vector by its own magnitude (to produce a unit vector) and multiply that by the force magnitude. I suspect that trigonometry is actually more efficient because you can use a look-up table.

EDIT: And what TythosEternal said.

Edited August 20, 2014 by Mr Shifty

[linkxsc]

Well from my terrible videogaem coding experience. (My end year project for the class was a "asteroids" type game, where you were in orbit. You had to keep the asteroids from hitting the planet or lose points. If you tried hard though you could use a tractor beam and redirect the roids into orbit themselves and tell them to get mined for extra points, and later turn them into little defense platforms that would help shoot down roids (they'd wait and only fire at where your mouse was pointed). All in all it was a fun little derp made for a project.

In that I did have n-body gravity dragging you around so you could get yourself in a semi orbit around an asteroid and such. I had planned to do more with it, but then I got a job and all that stuff ended.

Now the n-body worked just fine because there were basically no graphics or anything to work with. But as you started shooting stuff up and the asteroid count climbed, well on the school computers it started getting laggy around ~70 objects (cause remember, thats 70 objects interacting with eachother every frame, and lots and LOTS of math going down) After that I patched out the N body and just made them all check gravity to the planet, so then you could easily have hundreds of asteroids about without any lag. (ofcourse I'm sure there were countless things I could have done to streamline it more.) Game didn't bother at any point trying to track your own orbital path though. Just kinda had to speed stuff up until you were happy and then hope they didn't crash.

[Kermunist]

... or you can divide the relative position vector by its own magnitude (to produce a unit vector) and multiply that by the force magnitude.

I have a feeling that we basically agree how this works, it's just I wasn't quite clear originally that my division when calculating the force was a vector division by scalar. We're both thinking of using the same operations, but I had considered the hard ones as part of the calculation of the force, and you considered them part of the vector addition.

I edited my last post before seeing your latest, it might be a bit clearer now anyway. 3 divisions and a square root (which I had forgotten about) are the things that will take the most time. Still quicker than you will get with trig I think.

[K^2]

2) I'm talk about a per-body cache. May not be necessary, depends on the complexity of the gravity calculation. I'll assume I'm missing some important fact about gravity.

Computing gravitational pull isn't the hard part. Caching it is overkill, and is likely only to make things more complex than they are, since you'll still need to compute corrections between the grid points.

I'm guessing you've never studied the topic of integrating equations of motion. Simply computing forces isn't sufficient. Forces are going to change from one frame to another, and if you do updates as if they are just spikes of impulse at each point (Forward Euler), you'll end up with huge error accumulated. This is true even without involving gravity. For that reason, video games typically use Velocity Verlet. But for gravity, even that is absolutely terrible. The gravitational potential is 1/r, which means no power expansion is ever going to be exact. Worse yet, in n-body problem, the bodies are in motion. Which means that there is no coordinate system (Such as barycenter coordinates for 2-body problem) in which energy is conserved. So even if you could construct a proper conservative integration technique, it would be useless.

What does that mean in terms of the game? That results of integration aren't going to converge. No matter how small time steps you take, over large enough time frame, the error will be too high. That wouldn't be too bad if you only had to do this once. I mean, it's a game, it doesn't have to have NASA-grade precision. But then any tiny changes in initial condition would result in different integration results. This could still be ok if no forces were ever applied to ships, but add any kind of force on the craft, and everything flies to hell. Did you notice that sometimes, when you have a very iffy intercept, the predicted trajectory starts to flicker between two states? Well, that's how absolutely everything in the game would be. Every single trajectory would jump all over the place, and maneuvers as simple as Mun intercept would be come almost impossible.

It's possible to improve on that. You can use very complex implicit integration techniques, similar to these that NASA uses, to greatly improve precision. This is where CPU intensity comes in. These are very heavy computations. And they can still only tell you where a particular asteroid or a comet will pass within a huge margin of error. That's ok for real space flight, and for hard-core simulators, because you never expect to just do a burn at LEO, and do a perfect insertion into Moon's orbit. You'll have to do additional burns along the way, as you see that what you predicted in original computation isn't matching the real world. A long interplanetary mission would have many scheduled correction burns which would be adjusted as tracking data is updated. But this would fundamentally change KSP's gameplay. It would take mission planning from something that everyone can do, to something where you'd need a master's in celestial mechanics to understand. Making the game too complex for its own good both in terms of internal code and in terms of gameplay.

There is a sort-of compromise that might be workable. Instead of true n-body problem, you can do a perturbed central potential problem. This would involve computing dynamics of the orbital elements that describe the conic section instead of dynamics of the actual craft. This would give you a lot of interesting features. For example, it'd be entirely possible to have Lissajousand Halo orbits around Lagrange points. It'd also allow for things like orbital precession due to the oblateness of the planets. A lot of interesting, fun things, which would still allow for fundamental gameplay to be pretty much the same.

It is still way, way more complex computationally. And I'm not sure about stability near edges of Hill spheres. Could make for an interesting mod, however. I'll have to look into it.

[linkxsc]

Also dont forget that the game would have to be doing all these gravitational computations on top of the already rather straining (to some machines) inter-part forces within a craft itself. (this might be passed over perhaps if there were a mod that after "building" your ship locked it in as a rigid structure, but then thats kinda ruining the fun of the game)

In an extreme case, adding in full n-body physics would also being to cause tidal forces within your ship structure. (ship is a certain semi rigid shape moving at a constant speed. The parts of the ship closer the planet are moving faster when referenced from the planet than the parts farther out, and this could induce forces in the ship, at least I think it would. Don't really remember the rules for that, but I'm pretty sure the net effect would be your ship spinning wonkily, although on the scale of the game and by extension IRL it's not too much of a problem for anything we've built yet).

[andrewas]

IIRC, they took tidal forces out of the game, gravity is calculated once for the CoM of the ship rather than calculated independently for every part. For the most part, this has no observable effect on anything, it does mean you can't do gravity gradient stabilization but since ships on rails don't rotate that's pointless anyway.

[K^2]

Linkxsc, you can do rigid body and still do stress calculations. That would already be faster, provide more stable flight, result in fewer errors, and be more aesthetically pleasing. It would not reduce the design challenge, and unless you are in it for the absurd bouncing contraptions, I don't see how it would take away from the fun of it either. That is another mod idea, by the way, but it would require some highjacking, since this would not be covered by general plugin stuff.