Previous ToC Up Next

3. The NBody Class

3.1. Code Listing

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

3.2. Starting Up

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-3
An N-body system is read in, and all that happens is a single call to evolve, while the command line arguments are being passed along, just as we already saw in the previous code.

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.

3.3. Finding the Next Particle

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.

3.4. Output

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.

3.5. Wordlines and Worldpoints

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.

3.6. Verification

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 = 2
We 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-06
A 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-13
Okay, 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-12
Perfect, 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
Bob: There are larger than I would have expected. Especially the constant time step run is off by an amount that is two orders of magnitude larger than the relative energy conservation.

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-08
And 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
Alice: Now they are indeed all comparable. Good! I am beginning to believe that we may have done the right thing.
Previous ToC Up Next