Constant Acceleration Population

I couldn't sleep this morning. Woke up around 4am or so and just couldn't get my mind back to rest. That will make for a less than ideal day, but it had one significant positive. My brain was working on some useful ideas. The main one, and the one that finally pulled me out of bed to code it up, was to add a constant acceleration to my BasicPopulation. That sentence probably needs a bit of explanation.

My simulation code is rather modular. One of the more significant modules is an element called a population. The population keeps track of the bodies in the simulation and has logic in it to handle their movement, find collision times, process collisions, etc. I basically have two populations that I use. My ring simulations use what I call the GCPopulation. GC stands for guiding center. In this population, the particles don't move on straight lines. Instead, they follow guiding center motion which is a first order analytic solution to Hill's equations. This complicates collision finding, but has benefits that make it well worth while otherwise. Basically, the gravity of Saturn is implicit.

For general research purposes when I want to look at other types of systems I have created a "BasicPopulation". In this population particles travel on straight lines between collisions in the middle of a time step. I was using this for the simulations for Josh and I had a second force of constant acceleration in addition to the collision forcing.

This constant acceleration causes some issues numerically that are best resolved by taking small time steps. This is common in numerical integrations. The integration is more accurate if you take smaller time steps. The simulation also really grinds at the end when the particles are sitting right on top of one another in a close packed configuration. Part of the reason it grinds in this configuration is that the gravity force is applied as a kick at the beginning of the time step and all the particles are given a downward motion which then has to be canceled out by collisions propagating back up through them.

The thing about constant acceleration is that you really don't need a numerical integrator. There is a nice analytic solution to the motion, it is a parabola. So the idea was to create a population that uses that motion somewhat like the GCPopulation uses guiding center motion. When I turned this around in my head a bit (and didn't fall back asleep as a result) I realized that I didn't even need a new population. It's an easy addition to the BasicPopulation. How easy? Turns out to be one line. It was a little harder to make it so one other element of the simulation would work, but that wasn't too bad either. Now I'm testing it.

Why do I care? My hope is that at this will have a very helpful impact on the end of the simulation when all the particles are tight packed. Without the constant acceleration force the particles won't all be given a kick in the downward direction. Instead, they will follow a true parabola with a derivative of zero at apex. This should provide much nicer numerical behavior at that point and allow me to take larger time steps because I don't have to worry about numerical accuracy of the integration.

Now if I can just get that locking graph completed so I can drop a marble on the final configuration.

Two Argument Recursion

So today Cameron Swords and I talked through a possible solution to the problem of the last post. The solution involves using the tree to build a graph of "lock nodes". These are nodes in the tree that we use as locks instead of having a grid of locks to tell us if a collision is safe. In theory this method should work quite well. It will be as dynamic and optimized as the tree and have little additional overhead. For simulations with a size distribution of a dynamic distribution of particle relative velocities, it should provide optimal locking to keep the simulation running nicely on multiple threads.

One little detail of this method that I felt was worth commenting on was the creation of this graph. The "lock nodes" need to have some basic properties. There should be one on the path between the root and any leaf. It could be the leaf. (In theory it could be the root, but that would suck because there would be no parallelism in the collisions.) We would also like to use nodes that are as far down the tree as possible without making the graph too dense. The edges in the graph connect nodes that are sufficiently close that it isn't safe to do collisions in both at the same time. So we want nodes that have a size on the order of the search distance for collision pairs. Identifying these nodes isn't that hard. It can be done with a nice recursive decent. Once we have them we would have to build the graph which could be done easily in code by doing a recursive decent for each of the lock nodes. Such an approach is easy to describe and would be O(n log n) with a coefficient that isn't huge so it wouldn't be the limiting factor in the simulation. However, there is a better way.

This better way, as the title of this post alludes to, is to do recursion over two arguments. In this way, we can identify the lock nodes and build the graph all in a single recursive call. However, that recursive call now recurses on two arguments and basically runs through pairs of nodes. The advantage is that it prevents it from doing work on the top potion of the tree multiple times and should provide something close to that absolute minimum work needed to identify all of the pairings needed for the graph.

This isn't the first time I have used this type of method. In fact, it is the fourth time and the third time in this particular simulation code. I first used the method for my implementation of a symplectic tree code that I have worked on with Hal Levison. In that code the interactions occur only between nodes on the same level of the tree. As such, the two argument recursion makes great sense. I had never seen it used before and still haven't, so for all I know, it hasn't been widely discovered as an approach for numerical work.

My second use was in the collision finding in my main simulation code when a tree is used for finding the collisions. More recently I updated gravity to use it as well. In the case of gravity, switching to this approach dropped the number of distance calculations and improved code performance by a factor of 2-3 times. Now it looks like I'll be using it a fourth time to find pairs of lock nodes that should have edges between them in the gravity/collision tree.

As a general method, what this approach does is to efficiently run through and identify pairs of nodes in a tree when there is a criteria that can eliminate subtrees from consideration. In my work, the criteria is a distance issue, but it could be other things. This could also be generalized to recursion over three or more arguments for finding triples, etc. What makes it efficient is that the recursion doesn't pass over the top level nodes doing the check for inclusion multiple times. Instead, it goes through them once for each argument that it is recursing on and then goes on down to other pairs where the check is more significant.

(Interested readers might wonder how one could make a single recursive call like this operate in parallel. If anyone requests a description of how we have found to do that in a way that leads to generally good load balancing, I'd be happy to post it at that time.)