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.
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:
applyTreeForceshas 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?
applyTreeForces()return the maximum
DvDtitself, 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.
dtDynmethod in an unreachable if branch (lines 404-427
TreeGravity.cc). This can be cleaned up when we feel it's safe.