Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo


#4 branch grav_dt_mods ready for review and merge

Naor Movshovitz

Branch ready for review and merge. In this branch I modified the time step vote of the TreeGravity object, as follows.

Concept. I implement the dt vote based on the most strongly accelerated node. Specifically, if the maximum acceleration experienced by any node is DvDt_max, then I take the time scale to be t_scale = sqrt(system's bounding box length / DvDt_max), and the time step to be some specified fraction of that, using the settable mftimestep. This is the only time scale I use. The previously used dtDyn = 1/sqrt(G rho) is in fact included here. If, as is the case for a star/planet in equilibrium, the most accelerated particle is on the surface, its acceleration would be ~M^2/R, and the bounding box length would be 2R. This results in just that same time scale as 1/sqrt(2*G rho), if rho is interpreted correctly as the system's average density.

The pairwise r_relative/v_relative time constraint is removed. I thought at first this was an important constraint, but now I believe it is not. If two nodes want to swoosh past each other, let them. Whatever their relative velocities, if they are not strongly accelerating there is no reason to limit their time step.

Implementation. The way this can be done with minimal changes to source is to calculate the minimum dt, as before, inside TreeGravity::applyTreeForces(). But the dt is now updated in the outer loop, over node_i. First each node's DvDti is fully calculated. Then that DvDti is used to furnish a dt which is checked against the current dt to select a minimum.

For example, with my polytropic pseudo Jupiter, the maximum DvDti is ~23 m/s^2 - indeed the surface gravity of jupiter. And the bounding box length is ~1.4e8 - the diamater of jupiter. Resulting in a time scale of ~2500 s, indeed just what 1/sqrt(2*G*rho) would suggest. But this method would yield a good time step even when the planet is severely deformed.

Some details we to review before merging:

  • Did I miss some parallel aspect of this new calculation? applyTreeForces has some parallel structure that I don't fully appreciate, and presumably this is why we don't need a call to allReduce to get the global mindt?
  • If refactoring is safe and easy, I think it would be more helpful to have applyTreeForces() return the maximum DvDt itself, and have TreeGravity::dt() use that to derive a time step, just once. This would work since the length scale is now the global mBoxLength, and maybe save some cpu cycles.
  • I left the old dtDyn method in an unreachable if branch (lines 404-427 TreeGravity.cc). This can be cleaned up when we feel it's safe.