|
Orrery
A GPU-accelerated N-body gravitational simulator
|
Everything the simulator does, reachable from Python, with the particle state as NumPy arrays that share memory with the run rather than copies of it.
From a clone:
That builds the whole C++ library, which takes a couple of minutes, and installs the orrery package with the extension inside it. The only requirements are CMake 3.25 or later and a C++20 compiler; pip fetches the rest into an isolated build environment.
For the example notebooks:
A wheel is the wrong tool for that, since every edit would mean a reinstall. The python preset writes the extension into a directory that is already a working package:
build/python/python then holds an orrery directory with the extension and the package sources in it, so putting that on PYTHONPATH gives a working import against the build tree:
The test suite is pytest and runs under CTest with everything else, so a single ctest invocation is the whole definition of done rather than most of it. It needs NumPy and pytest; a configure that cannot find them says so and leaves the suite unregistered rather than failing it.
Configuration is the same record the configuration file parses into, so a run set up here and a run set up from a document are the same run. write_configuration turns one into the other:
and orrery run cluster.orrery then does what the script did.
This is what the bindings are for. The component arrays are NumPy views of the storage the solver is writing into:
Positions are held as three contiguous arrays rather than as one array of triples, because that is what makes the force kernel fast (ADR-0004), so an (N, 3) array of them exists nowhere in memory. The interface says so:
| What you want | How to get it | Cost |
|---|---|---|
| One component | particles.position_x | A view. No copy |
| All three | orrery.components(particles) | Three views in a tuple. No copy |
| An (N, 3) array | orrery.stacked(particles) | A copy |
| Distances from the origin | orrery.radii(particles) | A copy of one array |
The ten arrays are position_x, position_y, position_z, velocity_x, velocity_y, velocity_z, acceleration_x, acceleration_y, acceleration_z and mass. ADR-0040 records why there is no (N, 3) property.
A ParticleData you built is writable, and assigning into its arrays sets up a configuration:
A running simulation's state is not. simulation.particles is a ParticleView whose arrays have NumPy's writeable flag cleared, because the integrators require the acceleration array to hold the acceleration at the current positions on entry to every step and writing into a live run would break that silently. To change a run's state, build the state you want and hand it over:
restore re-establishes the invariant, at the cost of one force evaluation. ADR-0041 records the decision.
A view keeps the object that owns its storage alive, so an array outlives the expression that produced it and outlives the last Python reference to the store. It does not survive a reallocation: resize, reserve, add and Simulation.restore may move the component arrays, and a view taken before one of those points at freed memory afterwards, exactly as a std::span would. Take the view again after any of them.
compute_accelerations evaluates the forces on a set of particles once, in place, with whichever solver the configuration names, and returns the work it took. It is how an approximation is measured against the reference:
The interaction counter reports what the algorithm did rather than how long the machine took, which is what makes that comparison a statement about the methods.
Simulation.step, Simulation.run, Simulation.measure and compute_accelerations release the interpreter lock for the duration of the work, so the solver's own threads have the machine to themselves and another Python thread can watch a run while it happens.
assemble releases it for the sampling and the first force evaluation, which is where its time goes.
orrery.dtype is the NumPy scalar type the library was built with, float64 by default and float32 in a single-precision build. A tolerance in a script should be written in terms of it rather than assumed, since one written for double precision silently asserts nothing when the same code runs against a single-precision build.
In python/notebooks/. They are committed without outputs, so running one is the only way to see its figures, and every claim they make is asserted before it is plotted.
| Notebook | What it shows |
|---|---|
| 01_validation.ipynb | The measured convergence order of each integrator, bounded against secular energy error over four hundred orbits, and a sampled Plummer sphere against the closed-form model |
| 02_solver_accuracy.ipynb | Barnes-Hut against direct summation: error against opening angle, error against cost, the interaction saving as N grows, and the momentum conservation a tree gives up |
| 03_galaxy_collision.ipynb | The demonstration scenario, run and plotted, with the energy it conserved and where the material ended up |
Run them from a clean environment with:
or execute them without opening anything:
The numbers they produce are the numbers the C++ test suite asserts. The measured convergence orders come out at 1.9998, 4.0006 and 4.1659 against stated orders of 2, 4 and 4, which are the figures in the README, taken by a different route.
The trajectory and checkpoint readers, and the renderer. Neither is an oversight.
A run driven from Python reports through the NumPy views of its own state, which is strictly more than a binary file it would have to parse afterwards, so assemble attaches no output and the output section of a configuration is ignored. The command-line program is what writes trajectories and checkpoints, and the formats are specified in docs/formats/ well enough to be read by something other than this program.
The renderer needs a window, an OpenGL context and a display, none of which belongs in an extension module. orrery-view is what watches a run.