How can I reduce the amount of storage needed for a gravitational n-body simulation?

Question

I am currently attempting to create a gravitational n-body simulation using a modified Barnes-Hut algorithm, to be more amicable to GPU computation. This is primarily as a learning project. My goal is to simulate a number of stars comparable to that in real galaxies, meaning on the order of hundreds of billions to tens of trillions, but even a few million would be useful. It is very unlikely that I will be able to compute this at a speed amicable to display, meaning that I must pre-compute the data and look at it after the computation finishes. To do this, my first idea for how to store the data is to create a file that has the locations of all the stars concatenated together for each moment of discretized time, and then the next moment concatenated to that, to make something like the following, where the data in each bracket represents a single frame:

{x₁, y₁, z₁, x₂, y₂, z₂, …. xₙ, yₙ, zₙ}, {x₁, y₁, z₁, x₂, y₂, z₂, …. xₙ, yₙ, zₙ}, {x₁, y₁, z₁, x₂, y₂, z₂, …. xₙ, yₙ, zₙ}

Alternatively, this format can be described in C++ pusedo-code (portability between C++ implementations is not important):

void writeData(std::vector<Frame> frames, std::ostream &out){
    //decoding knows how many points there are in each frame as a property of the file format, so it    can read an entire frame at a time until EOF
    for(const Frame &frame : frames){
        for(Point &point : frame.points()){
            float x = point.x();
            float y = point.y();
            float z = point.z();
            out.write(&x, sizeof(float));
            out.write(&y, sizeof(float));
            out.write(&z, sizeof(float));
        }
    }
}

The size of the data given by this format in bytes is n*3*2*60*25 for n particles, one minute of data, 25 frames per second, and half precision floats. For one billion particles, this works out to 16 terabytes for one minute of video with one billion particles, something that I physically don’t have even close to enough hard drives to store (and I have a lot of hard drives). I also doubt that this data will compress well with standard lossless compression algorithms such as zlib, as it is fairly unstructured from a binary perspective. I also can’t think of any reasonably simple compression algorithms that would work well for this data.

Of course, I can encode a rendered frame of the particles in a 2-dimensional image and create a video with any of the many modern compression algorithms, but this sacrifices the ability to change the camera’s location while observing the data, something essential to gaining a good understanding of the three dimensional layout of the simulated particles (ie. consider that constellations appear to be completely different as Earth moves through the galaxy, over geologic time scales). Encoding multiple videos from different perspectives somewhat mitigates this, but nothing compares to the ability to control the camera’s location and angle during playback.

How can I enable camera movement with pre-computed data and a large number of particles without having thousands of dollars to drop on massive hard drives? I think that this is possible because regular, two dimensional, video is a similar problem that is now solved enough to be practically useful, and because other people have done n-body simulations (as I have seen videos of what can only be n-body simulations in news publications).

Answer 1

Float Warning

Be aware of floating point precision. The mantessa is only x bits long, which means that for values close to zero they can express very small distances, but at galactic distances the error can be quite significant.

It also makes a number of logical compression techniques harder.

Fixed Point Deltas

If you were to go with fixed-point numbers you could store deltas. This is based on the observation that most objects will not have traveled very far from the last frame. So instead of storing 3*N bytes per co-ordinate you could get away with say 3 or 6 bytes.

This will complicate your encoding scheme a little as you'll need to signal in advance the length of the deltas for the next k co-ordinates. But it should give a significant pay-off in terms of size.

Error-Corrected Guessing

Given that these points are moving, you can "guess" where the next point will occur. So the first two or three co-ordinates are full length, or full-length plus offset co-ordinates. The following co-ordinates after that are guessed using a line, or curve that fits the previous X points. The co-ordinate information written down is the error between that guess, and the actual co-ordinates generated by your simulation.

Frame of reference

Given that the objects are moving in a direction, you could reframe the deltas from global co-ordinates to local co-ordinates. In local co-ordinates the same or similar deltas are likely to occur improving compression efficiency generally.

Similarly in local co-ordinates the previous delta can be considered a form of guess about the next location, with the delta correcting the guess and consequently being smaller.

Encode a physics delta, such as acceleration

Considering that you only seem to be worried about location, you could take advantage of the naturally built in deltas that physics offers, such as acceleration and velocity.

Newtonian location/velocity/acceleration is relatively simple to calculate on the consuming end (and is sufficient for producing a movie). The simulation can workout the actual problem using whatever physics description and describe this on a global co-ordinate system. Then this is reverse described using Newtonian calculations to produce appropriate deltas for storage and later rendering.

In this sense the file starts with a location + velocity. Every subsequent co-ordinate is then derived from adding an acceleration to that dynamic and updating it for the next frame. The acceleration can be encoded raw, or a delta, or an error-correction to a guess - whatever works well.

Granularity

Generally many of objects will be moving along similar paths, you could find the mean/mode/avg movement for a group of co-located objects and describe that operation (translation + rotation) at a high-level, then describe the error between the global/group level "guess" at the new co-ordinates and the actual co-ordinates.

If the movements are consistent enough this will generally give much smaller deltas.

Unary and bit packing

If you manage to squeeze the deltas above down to frequently very small deltas/error corrections you can make use of unary and bit-packing to shrink it even smaller.

Unary is a number encoding not to dissimilar to holding up fingers to count.

0           //1 bit long, represents the number 0
10          //2 bits long, represents the number 1
11110       //5 bits long, represents the number 4
11111111110 //11 bits long, represents the number 10

You can use these unary numbers either directly to represent tiny deltas, or as the exponent to a formula that describes how many bits go into the co-ordinate deltas, etc... ie. 2*2^k would indicate co-ordinates that are: 2bits, 4bits, 6bits, etc...

By bit-packing in a 64-bit word you can squeeze in 10 - 2bit co-ordinate pairs.

0123456789ABCDEF 0123456789ABCDEF 0123456789ABCDEF 0123456789ABCDEF
110-Q1Q2Q3R1R2R3 S1S2S3T1T2T3U1U2 U3V1V2V3W1W2W3X1 X2X3Y1Y2Y3Z1Z2Z3

//the first 3bits is unary describing that the 64-bit record is 2bit co-ordinate triples.
//the - is an unused bit.
//Q1Q2Q3 through Z1Z2Z3 are the co-ordinate triples.

Mileage will vary though, so analyse some of the data your simulation produces and how similar it is between time-steps.