Previous | ToC | Up | Next |

** Alice**: I think I understand now how the ` Body` class works. Could you
also print out the whole ` NBody` class?

** Bob**: Sure. It's not very long, about the same size as the ` Body` class.
I keep being surprised at how compact Ruby code is.

class NBody def startup(dt_param) @e0 = ekin + epot @body.each do |b| b.pred_pos = b.pos b.pred_vel = b.vel end @body.each do |b| b.acc, b.jerk = b.get_acc_and_jerk(@body) b.time = @time b.next_time = @time + b.collision_time_scale(@body) * dt_param end end def evolve(c) @nsteps = 0 startup(c.dt_param) write_diagnostics t_dia = @time + c.dt_dia t_out = @time + c.dt_out t_end = @time + c.dt_end acs_write if c.init_out while @time < t_end np = find_next_particle @time = np.next_time if (@time < t_end) np.autonomous_step(@body, c.dt_param) @nsteps += 1 end if @time >= t_dia sync(t_dia, c.dt_param) @nsteps += @body.size write_diagnostics t_dia += c.dt_dia end if @time >= t_out sync(t_out, c.dt_param) # we are now syncing twice, if t_dia = t_out @nsteps += @body.size acs_write t_out += c.dt_out end end end def find_next_particle next_time = VERY_LARGE_NUMBER next_particle = nil @body.each do |b| if next_time > b.next_time next_time = b.next_time next_particle = b end end next_particle end def sync(t, dt_param) @body.each{|b| b.forced_step(@body, t, dt_param)} @time = t end def ekin # kinetic energy e = 0 @body.each{|b| e += b.ekin} e end def epot # potential energy e = 0 @body.each{|b| e += b.epot(@body)} e/2 # pairwise potentials were counted twice end def write_diagnostics etot = ekin + epot STDERR.print <<END at time t = #{sprintf("%g", @time)}, after #{@nsteps} steps : E_kin = #{sprintf("%.3g", ekin)} ,\ E_pot = #{sprintf("%.3g", epot)} ,\ E_tot = #{sprintf("%.3g", etot)} E_tot - E_init = #{sprintf("%.3g", etot - @e0)} (E_tot - E_init) / E_init = #{sprintf("%.3g", (etot - @e0)/@e0 )} END end end

** Alice**: We already talked about the generic case, with a call to
`autonomous_step`. Now I would like to see how the whole
thing starts up.

** Bob**: At the very bottom of the file, the ` evolve` method is evoked
as follows:

: inccode:.nbody_ind1.rb-3An N-body system is read in, and all that happens is a single call to

The startup function makes sure that the acceleration and jerk are being
computed properly at the beginning of the integration. But in order to
do so, it first has to set the predicted positions and velocities
`pred_pos` and `pred_vel` to the right values, namely the
values of the initial positions and velocities.

In the ` Body` class we saw the following generic order: 1) we predict
the positions of all particles, using `predict_step`; 2) we
compute the acceleration and jerk on one particle and set the next
time step value; 3) we step that particle forwards using `correct_step`.

We use the same order here, at startup, which occurs at the time
`@time`, the instance variable for the ` NBody` class. This is
the time that has been read in from the initial snapshot; and if there
is no time listed in that input file, the default time `@time = 0`
is used default, since this is the time that has been assinged by the
`NBody#initialize` method in file `nbody.rb`.

For startup, step 1), prediction, is trivial: starting at time `@time`
we predict where all particles will be at time `@time`, in other words
we don't have to do anything; we just copy the position and velocity values
that were just read in.

Step 2), the calculation of acceleration and jerk and the next time step
value, can only be done after step 1) has been completed for all particles.
This is the reason that we need to different ` each` loops in ` startup`.

Step 3), correcting the positions and velocities, is not necessary at startup: no particle has been moved, so there is nothing to correct.

** Alice**: Thanks, that's very clear. Then, after startup, we enter the
` while` loop in evolve, evoking `autonomous_step` for each next
particle that needs to be propagated. Let's have a look at how
`find_next_particle` is implemented.

def find_next_particle next_time = VERY_LARGE_NUMBER next_particle = nil @body.each do |b| if next_time > b.next_time next_time = b.next_time next_particle = b end end next_particle end

** Bob**: Here we loop over all particles, to find the particle with the
earliest next sales date `next_time`. The method
`find_next_particle` then returns the body that has the smallest
`next_time` value.

** Alice**: So time is marching forward. When we enter
`find_next_particle`, the ` NBODY` time `@time` has a certain
value. We then are given a body ` b` which has a value `b.next_time`
for which we know that `b.next_time > @time` and we also know that
for all other bodies ` ob` we have `ob.next_time >= b.next_time`.

** Bob**: Correct.

** Alice**: This means that in the interval between `@time` and
`b.next_time` there is nothing that needs to be done. Then we
set the ` NBody` time `@time` equal to `b.next_time`, and
we push particle ` b` forwards by one step.

Having done that, we again look for the next particle, through a new
call within ` evolve` to `find_next_particle`. Perhaps the time
for the newly found particle is the same as that for the previous
particle, in which case `@time` doesn't change; perhaps it is
later, in which case `@time` is updated to the latest sales
date of the new particle.

Okay, I think I got it! I just had to reconstruct the steps for myself. This is not a trivial algorithm!

** Bob**: As always, it is all pretty clear once you understand it clearly,
but I must admit, it took me quite a while to figure it out, the first
time I came across an individual time step code.

** Alice**: We're almost there. I just have to figure out now how the
particles are synchronized before you do an output, either a full
particle output in terms of a snapshot, at time `t_out`, or a
diagnostics output at time `t_dia`.

Let me put the need for synchronization in perspective. In the case
of our constant time step code `nbody_cst1.rb`, synchronization
was trivial: in typical usage, you specify output times that are multiples
of the prescribed time step. And if you were to set a time step to, say,
0.01 and an output time to an incommensurable time 1/3, then the code would
slightly overshoot, and give the output at time 0.34, instead of time
0.333333. We could have added an option to allow us to shrink the
last time step in such a case, to reach the value 1/3 exactly, at
least within double precision. However, there was no pressing reason
to do so.

In the case of our shared time step code `nbody_sh1.rb`, we did
build in such an option. Invoking it by default would let the code
overshoot, since in general the chance is virtually zero to get
commensurability between a prescribed output time and dynamically
adjusted time step values. However, by invoking the code with the
option `nbody_sh1.rb --exact_time`, we could force the code to
halt at the exact time desired. This was essential for us in order to
measure the phase space distances between different N-body systems
generated in the output of different runs.

Now in the case of our individual time step code `nbody_ind1.rb`,
we have no choice. We cannot afford the luxury of skipping synchronization,
since in that case we would get an output of a bunch of particles, all at
different times. We couldn't even compute the energy of sych a system!
And this is the reason that you synchronize, for every type of output.

** Bob**: Yes, that is the big picture. Note, however, that my insistence
on synchronization is not the only solution. I wanted to get
something coded up quickly, and I knew it was essential for measuring
energy conservation, so I implemented syncing by default. However, if
you just want to make a dump of the system, to allow future restarts,
there is no reason to insist on all the particles having their state
set at the same time.

** Alice**: Ah, yes, of course, I had forgotten that. Indeed, if you
* do* synchronize, you are guaranteed * not* to follow the same trajectory
when you restart from an earlier output. We showed this in the previous
volume, in the case of a shared time step code.

Well, what do you think, shouldn't we implement such a `ragged' output, with the orbits of different particles extending in time to different distances?

** Bob**: Ragged?

** Alice**: I'm trying to visualize what is going on. I can see a picture
in front of me, in spacetime, with each particle following a worldline,
where some worldlines are extended further than others, like a ragged
carpet where the different strands have different lengths.

** Bob**: You always manage to complicate things! As soon as you understand
how something works in three dimensions, you add another dimension for
good measure!

** Alice**: I'm serious, though. I do think it may be an interesting way
of looking at the evolution of an N-body system. Normally we
associate a spacetime view with special relativity, where particles
follow worldlines in space and time, but there is no reason not to use
such a picture in classical mechanics as well.

** Bob**: So instead of talking about the N-body problem, from now on you
want to talk about the N-worldline problem?

** Alice**: Well, why not?

** Bob**: I know why not. Nobody will understand what you're talking about.

** Alice**: We'll worry about that later. Why not pursue this idea for a
moment. Rather than thinking in terms of a bunch of separate bodies,
as isolated point masses in 3-D, it may be more natural to think of
them as strings in 4-D. Such a picture could actually be helpful in
visualizing the basic idea of an individual time step code.

** Bob**: A string theory for classical gravity?

** Alice**: If you like! But without Calabi-Yau manifolds, don't worry!
We'll try to keep it simple.

** Bob**: What are Calabi-Yau manifolds?

** Alice**: Something quite a bit more difficult to visualize than
individual time step algorithms. Forgot about that.

** Bob**: With pleasure; life is already complicated enough.
So instead of having N particles in space, you want us to look at N
worldlines in spacetime. And each worldline has a number of knots,
one at each point where we compute the positions and velocities.
Do you want to call those knots worldpoints?

** Alice**: That would be a natural term, although I haven't heard it before.
But why not? In special relativity textbooks, these points in 4-D are
called events, since they are associated with a particular time and place.
However, the word worldpoint is more descriptive, in that they convey the
notion that they belong to a particular worldline. Different worldpoints
on the same worldline then share the same `body_id` identity, but
are different in having different values for the time `@time`

** Bob**: Taking a time step is then a temporal operation, while syncing
for output is a spatial operation?

** Alice**: You're getting it! See how natural it is to view the N-body
problem in four dimensions?

** Bob**: The N-worldline problem.

** Alice**: I'm happy to keep calling it the N-body problem, not to worry.
And I do think this way of looking at things has mileage. Synchronization
means constructing a cross section of the full bundle of worldline with a
hypersurface at a constant time.

** Bob**: And what is the mileage in being abstract to that degree?

** Alice**: Let me think about it a bit more, and continue tomorrow.

** Bob**: Good idea: with all this talking about spacetime, we forgot the
time. It's getting late.

** Alice**: Ah, but before we forget, did you test your individual time step
scheme, to make sure it worked?

** Bob**: I did, but just to check, and to show you it's in good shape, let's
do a quick run.

** Alice**: Let's take all three codes we have now, for constant, shared,
and individual time steps. And let's give them a really good workout.
Shall we shoot for an accuracy of ? If we get too
close to machine accuracy, it may be more difficult to interpret our
results.

** Bob**: Fine with me. Let's start with the same initial conditions:

|gravity> kali mkplummer.rb -n 4 -s 1 | kali nbody_set_id.rb > test.in ==> Plummer's Model Builder <== Number of particles : N = 4 pseudorandom number seed given : 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 actual seed used : 1 ==> Takes an N-body system, and gives each body a unique ID <== value of @body_id for 1st body : n = 1 Floating point precision : precision = 16 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2We will have to try a few different time step accuracy parameters, to get roughly the right relative energy conservation. This may take a few repetitions. I'll start with constant time steps:

|gravity> kali nbody_cst1.rb -t 1 < test.in > tmpc.out ==> Constant Time Step Code <== Integration method : method = hermite Softening length : eps = 0.0 Time step size : dt = 0.001 Interval between diagnostics output: dt_dia = 1.0 Time interval between snapshot output: dt_out = 1.0 Duration of the integration : t = 1.0 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 at time t = 0, after 0 steps : E_kin = 0.25 , E_pot = -0.5 , E_tot = -0.25 E_tot - E_init = 0 (E_tot - E_init) / E_init = -0 at time t = 1, after 1000 steps : E_kin = 0.248 , E_pot = -0.613 , E_tot = -0.365 E_tot - E_init = -0.115 (E_tot - E_init) / E_init = 0.461 |gravity> kali nbody_cst1.rb -t 1 -c 0.0001 < test.in > tmpc.out ==> Constant Time Step Code <== Integration method : method = hermite Softening length : eps = 0.0 Time step size : dt = 0.0001 Interval between diagnostics output: dt_dia = 1.0 Time interval between snapshot output: dt_out = 1.0 Duration of the integration : t = 1.0 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 at time t = 0, after 0 steps : E_kin = 0.25 , E_pot = -0.5 , E_tot = -0.25 E_tot - E_init = 0 (E_tot - E_init) / E_init = -0 at time t = 1, after 10000 steps : E_kin = 0.313 , E_pot = -0.563 , E_tot = -0.25 E_tot - E_init = -6.62e-07 (E_tot - E_init) / E_init = 2.65e-06A bit better than needed, but this will do. Now for shared time steps:

|gravity> kali nbody_sh1.rb -t 1 --exact_time < test.in > tmps.out ==> Shared Time Step Code <== Integration method : method = hermite Parameter to determine time step size: dt_param = 0.01 Interval between diagnostics output: dt_dia = 1.0 Time interval between snapshot output: dt_out = 1.0 Duration of the integration : t = 1.0 Force all outputs to occur at the exact times Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 at time t = 0, after 0 steps : E_kin = 0.25 , E_pot = -0.5 , E_tot = -0.25 E_tot - E_init = 0 (E_tot - E_init) / E_init = -0 at time t = 1, after 1927 steps : E_kin = 0.313 , E_pot = -0.563 , E_tot = -0.25 E_tot - E_init = -2.33e-10 (E_tot - E_init) / E_init = 9.32e-10 |gravity> kali nbody_sh1.rb -c 0.003 -t 1 --exact_time < test.in > tmps.out ==> Shared Time Step Code <== Integration method : method = hermite Parameter to determine time step size: dt_param = 0.003 Interval between diagnostics output: dt_dia = 1.0 Time interval between snapshot output: dt_out = 1.0 Duration of the integration : t = 1.0 Force all outputs to occur at the exact times Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 at time t = 0, after 0 steps : E_kin = 0.25 , E_pot = -0.5 , E_tot = -0.25 E_tot - E_init = 0 (E_tot - E_init) / E_init = -0 at time t = 1, after 6424 steps : E_kin = 0.313 , E_pot = -0.563 , E_tot = -0.25 E_tot - E_init = -6.12e-13 (E_tot - E_init) / E_init = 2.45e-12 |gravity> kali nbody_sh1.rb -c 0.002 -t 1 --exact_time < test.in > tmps.out ==> Shared Time Step Code <== Integration method : method = hermite Parameter to determine time step size: dt_param = 0.002 Interval between diagnostics output: dt_dia = 1.0 Time interval between snapshot output: dt_out = 1.0 Duration of the integration : t = 1.0 Force all outputs to occur at the exact times Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 at time t = 0, after 0 steps : E_kin = 0.25 , E_pot = -0.5 , E_tot = -0.25 E_tot - E_init = 0 (E_tot - E_init) / E_init = -0 at time t = 1, after 9635 steps : E_kin = 0.313 , E_pot = -0.563 , E_tot = -0.25 E_tot - E_init = -1.24e-13 (E_tot - E_init) / E_init = 4.96e-13Okay, that one's below our threshold too. Now for individual time steps:

|gravity> kali nbody_ind1.rb -t 1 < test.in > tmpi.out ==> Individual Time Step Hermite Code <== Parameter to determine time step size: dt_param = 0.01 Interval between diagnostics output: dt_dia = 1.0 Time interval between snapshot output: dt_out = 1.0 Duration of the integration : t = 1.0 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 at time t = 0, after 0 steps : E_kin = 0.25 , E_pot = -0.5 , E_tot = -0.25 E_tot - E_init = 0 (E_tot - E_init) / E_init = -0 at time t = 1, after 4257 steps : E_kin = 0.313 , E_pot = -0.563 , E_tot = -0.25 E_tot - E_init = 7.24e-10 (E_tot - E_init) / E_init = -2.9e-09 |gravity> kali nbody_ind1.rb -c 0.002 -t 1 < test.in > tmpi.out ==> Individual Time Step Hermite Code <== Parameter to determine time step size: dt_param = 0.002 Interval between diagnostics output: dt_dia = 1.0 Time interval between snapshot output: dt_out = 1.0 Duration of the integration : t = 1.0 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 at time t = 0, after 0 steps : E_kin = 0.25 , E_pot = -0.5 , E_tot = -0.25 E_tot - E_init = 0 (E_tot - E_init) / E_init = -0 at time t = 1, after 21282 steps : E_kin = 0.313 , E_pot = -0.563 , E_tot = -0.25 E_tot - E_init = -9.77e-13 (E_tot - E_init) / E_init = 3.91e-12Perfect, just on the mark.

** Alice**: And now let's check whether the results are really similar:

|gravity> cat tmpc.out tmps.out | kali nbody_diff.rb ==> 6N-dimensional phase space distance between two N-body systems <== Floating point precision : precision = 2 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 6N-dim. phase space dist. for two 4-body systems: 6.8757897588114396e-06 |gravity> cat tmps.out tmpi.out | kali nbody_diff.rb ==> 6N-dimensional phase space distance between two N-body systems <== Floating point precision : precision = 2 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 6N-dim. phase space dist. for two 4-body systems: 3.5201108698382295e-11

** Alice**: Yes, that is a bit much. The good news is that the shared time step
code and the individual time step code give results that correspond to about
. That may not be so bad; there is no guarantee that
the phase space distance should give the same result as the energy error.

** Bob**: Let me run the constant time step code with a two times smaller
time step size. That should increase the accuracy by more than an order
of magnitude, bringing all three codes more in line:

|gravity> kali nbody_cst1.rb -t 1 -c 0.00005 < test.in > tmpc.out ==> Constant Time Step Code <== Integration method : method = hermite Softening length : eps = 0.0 Time step size : dt = 5.0e-05 Interval between diagnostics output: dt_dia = 1.0 Time interval between snapshot output: dt_out = 1.0 Duration of the integration : t = 1.0 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 at time t = 0, after 0 steps : E_kin = 0.25 , E_pot = -0.5 , E_tot = -0.25 E_tot - E_init = 0 (E_tot - E_init) / E_init = -0 at time t = 1, after 20000 steps : E_kin = 0.313 , E_pot = -0.563 , E_tot = -0.25 E_tot - E_init = -1.35e-08 (E_tot - E_init) / E_init = 5.4e-08And let's now compare all three pairs:

|gravity> cat tmpc.out tmps.out | kali nbody_diff.rb ==> 6N-dimensional phase space distance between two N-body systems <== Floating point precision : precision = 2 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 6N-dim. phase space dist. for two 4-body systems: 3.5758322798815124e-07 |gravity> cat tmpc.out tmpi.out | kali nbody_diff.rb ==> 6N-dimensional phase space distance between two N-body systems <== Floating point precision : precision = 2 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 6N-dim. phase space dist. for two 4-body systems: 3.5757721030515277e-07 |gravity> cat tmps.out tmpi.out | kali nbody_diff.rb ==> 6N-dimensional phase space distance between two N-body systems <== Floating point precision : precision = 2 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 16 Incremental indentation : add_indent = 2 6N-dim. phase space dist. for two 4-body systems: 3.5201108698382295e-11

Previous | ToC | Up | Next |