

Gravity simulation on a desktop computer 
The current implementations of Octree and KDTree are unbalanced in the sense that the spacedivision is done purely on a divideby2 basis. This means that one cell can end up with the majority of the bodies, leading to an asymmetric or unbalanced tree structure. It should be possible to produce balanced variants of each of these such that a roughly equal number of bodies is placed in each cell. This technique applies particularly cleanly to the KDTree and was the basis of the original work done by Joachim Stadel.
The balanced Octree and KDTree implementations could be made available in GravSim as further options and tested and compared against the unbalanced variants.
Note that the alternative Hilbert curve / treebased space divider (as used in the Ridlerstyle accelerators) is imnplicitly balanced because the tree construction is done bottomup rather than topdown.
GravSim version 1.11 includes a simple Kepler orbit analyser for binary stars in elliptical orbits. This works by checking pairs of stars via the collision timescale calculation and generating Kepler orbital elements if the eccentricity is less than 1.0 (i.e. a bound orbit). The calculation is done in the traditional way via the eccentric anomaly, reversecalculated from the mean anomaly by NewtonRaphson.
The next step would be to extend the Kepler regularisation to eccentricities greater than 1.0 (i.e. hyperbolic orbits). The logic for reversecalculating eccentric anomaly from mean anomaly is similar to the elliptic case, but uses hyperbolic sines and cosines instead. However, it is more difficult to control the numerical errors, particularly with high eccentricities and close encounters.
A further case to include in the regularisation would be eccentricities close to 1.0 (i.e. parabolic orbits). These arise quite rarely in practice, but can be solved with Barker's equation in conjunction with the other orbital elements.
A final case to include in the basic Kepler regularisation scheme is the simple linear orbit. These can arise in practice as the first step in coldcollapse simulations (i.e. where all initial velocities are zero).
The Pairwise Interactor (as used in conjunction with treebased space dividers to produce Dehnenstyle accelerators) has proved problematic to parallelise in attempts so far. The reason is that OpenMP does not support nested parallel execution, which forces the code to make a single decision to parallelise at an appropriate point in the tree. The trouble is that the Octree and KDTree implementations as they stand are unbalanced, which means that typically one processor gets to do all the work while the others finish quickly and sit idle. What is needed is either balanced implementations of Octree and KDTree or a way of measuring the degree of imbalance in the existing implementations. Once that is done, the details of the parallelisation will be similar to the existing BruteForce pairwiseparallel algorithm.
The current implementation of the GravSim PredictorCorrector integrators is ported from ACS and uses the simple PEC mode. This does a prediction prior to gravity evaluation and then a correction afterwards. The use of the correction typically improves energy errors by roughly a factor of 10 (according to Aarseth), without having to do any additional gravity evaluations.
However, it stands to reason that having done the correction, we now know that the original gravity evaluation was taken at the incorrect location and so can be improved. Therefore there is an argument for doing the gravity evaluation (and hence the correction) again. The resulting sequence is predictevaluatecorrectevaluatecorrect or P(EC)^{2}. Theoretically you could keep repeating this process, but experience and the law of diminishing returns suggests that twice is sufficient.
Once implemented, the P(EC)^{2} mode will be doing 2 gravity evaluations per timestep rather than just one. The crucial test of its effectiveness will be to compare the results against the original PEC mode but with twice as many timesteps. On the one hand, the accuracy of the integrator will improve with the number of timesteps according to its order. On the other hand, numerical errors due to repeated closeencounters can be expected to increase as a result.
The recent widespread availability of graphics cards with general processing capability makes them an obvious target for future GravSim development. This represents an alternative approach to the use of tree codes  i.e. throw hardware at the problem rather than trying to reduce the scale of the problem itself. Feedback from other researchers suggests that a factor of 10 improvement may be possible with CUDAbased nVidia graphics cards.
It is envisaged that the resulting system could mixandmatch GPU bruteforce vs a CPU treecode on a pertimestep basis, depending on the number of bodies involved. Clearly Shared or Constant timesteps (which involve the maximum number of bodies in each timestep) would run well on the GPU. However, recent experiments with large models (100,000 bodies or more) suggest that even individual timesteps would benefit. In particular, the initial gravity calculation at the start of the simulation and the final timestep prior to each logpoint always involve the maximum number of bodies. Even some of the intermediate steps involve sufficiently many bodies to make the treecodes take a noticeable amount of time, which suggests that the GPU will be a useful alternative here too.
Depending on progress, it may eventually be possible to run the tree codes on the GPU as well, opening up the way for even larger models.
At some point in the future, it may make sense to integrate GravSim with specialpurpose hardware such as the GRAPE6. It is anticipated that the impact on the GravSim system design would be similar to the generalpurpose GPUs, although it is not clear that running treecodes on specialpurpose hardware is a realistic possibility.
The Nbody literature makes extensive references to Kustaanheimo and Stiefel (KS) regularisation as a better way of keeping a lid on numerical errors via the the use of a generalised fictitious timescale. This is a clear candidate for inclusion in GravSim via a new plugin interface that gives a choice of orbit analyser, retaining the simple Kepler approach as an option.
The Nbody industrystandard reference is "Gravitational NBody Simulations, Tools and Algorithms" by Sverre J. Aarseth, Cambridge University Press, 2003. This book covers many complex subject areas but the extended regularisation could be included in GravSim at some point in the future. That would take it beyond a simple collisional simulator that copes with binary stars to the full gamut of triple stars, chain systems and hierarchical systems. The quest for numerical accuracy and dealing with the pathological cases that would otherwise cause excessive errors is the main theme of this work.
GravSim already employs a very basic neighbours scheme to determine the next timestep. For numericonly integration, it calculates the minimum collision timescale on a pairwise basis and then groups bodies with similar timescales (within a factor of sqrt(2)) into a common timestep. For analytic integration with the Kepler orbit analyser, it calculates the 2 least collision timescales on a pairwise basis. If the minimum collision timescale represents an elliptical orbit, then it uses the next one for timescale determination (assuming unperturbed motion in the meantime).
ACNS is a generalisation of this approach to neighbour lists greater than 2. The idea is that the nearerneighbour interactions are evaluated more frequently than those further away, which leads to an overall improvement of around a factor of 2 according to Aarseth.
This is a possibiity for future inclusion in GravSim, in conjunction with hierarchical and chain regularisation. However, it is not clear at this stage whether ACNS is an alternative to the collision timescale calculation or whether it can add anything via combination.
The GravSim tree codes are currently based on a fixedorder cartesian multipole expansion, up to the quadrupole term. An alternative would be to include a multipole representation in spherical coordinates, which is then more easily generalised to an arbitrary number of terms. A tradeoff is expected between the increased accuracy and the additional processing required. Bearing in mind that the whole point of a tree code is to go faster than bruteforce, current experience suggests that the next 2 terms in the sequence (octupole and hexedecapole) should be beneficial, but after that we should expect diminishing returns. Note that a plugin interface would be the ideal solution here  retaining a choice between the spherical and cartesian coordinates.
The accuracy of the cartesian multipole expansion currently used in the GravSim tree codes can also be improved by incorporating the next 2 terms in the sequence in cartesian form. Note that this will affect the optimum setting of the minimum batch size parameter, given the additional computation required during tree construction.
Once the spherical multipole moments are available in GravSim, it would then be possible to implement (in conjunction with the Pairwise Interactor) the local multipole expansions that are characteristic of the FMM. Note that this contrasts with the Taylorseries local expansions advocated by Walter Dehnen and currently used in the GravSim Pairwise Interactor.
One suggestion arising from the June 11th talk at MICA Simulations is the use of the "Aarseth Criterion" as an alternative to the collision timescale calculation used in ACS and GravSim. It may be possible to engineer a plugin interface that gives a choice of the method used in this area (i.e. how to choose the timscale for the next step). Note that GravSim has recently made further use of the collision timescale concept with the Kepler orbit analyser for binary stars. The inclusion of the Aarseth criterion will likely involve significant changes to several parts of this code.
The speed of light is not currently modeled in GravSim. Gravity is immediately felt everywhere at the same time and so subtle effects such as the precession of the orbit of Mercury will not be seen. Inclusion of GR would have a significant effect on all calculations, yet it would only really make a difference in certain specific situations. The ideal solution would therefore be to have it available as an option.
GravSim is collisional in the sense that it uses pure unsoftened Newtonian gravity by default. However, it represents bodies as point masses that never actually collide. Bodies typically fly off in a different direction (hyperbolic orbit) due to a close gravitational encounter, or slip past each other quietly if softening is used.
Modelling of star collisions is very different to (say) rock collisions and would need separate treatment to do it realistically. This is a very complex subject area but is a candidate for future inclusion in GravSim.
Although GravSim can be used to model stars in clusters and galaxies, it does not attempt to model the stars themselves. In reality a main sequence star will eventually evolve into a red giant and the transition could affect the results of the simulation, particularly in terms of collisions. Since Stellar Evolution is available via thirdparty packages in quite a mature form, the most sensible approach would be to link with these rather than including SE directly in the GravSim code. The relatively loose coupling between Nbody dynamics and SE makes this attractive in performance terms as well.
There are at least 2 different ways of integrating GravSim with SE:
The GravSim predefined model files were generated with the Plummer code as ported from ACS. There is a range of other models defined in the literature, each of which has its own merits and could be incorporated into GravSim. This could extend GravSim's reach beyond simple clusters to elliptical and spiral galaxies.
Existing packages such as Nemo and Starlab already include a good selection of models. Supporting their file formats would be a sensible next step.
GravSim is currently quite limited in the areas of model and orbit visualisation. There are opportunities to improve the GLUTbased GravView viewer. However, a better approach may be to link with other visualisation packages and projects.
One interesting such project is the OpenSim viewer produced for the MetaInstitute for Computational Astrophysics (MICA) which has been setup in the virtual world Second Life.
The realism of a gravityonly simulation becomes questionable as soon as you get down to the level of dust and gas. Basically you need more physics beyond just gravity when these small bodies collide and electromagnetism dominates.
Stateoftheart cosmological simulations tend to include gas as well as stars, often modeled via a technique known as SPH. One of the leading packages in this area is GADGET2 developed by Volker Springel.
There are a variety of other ways to model gravity and one that is often used in cosmological simulation is PM, useful for extremely large numbers of spreadout particles.