Previous ToC Up Next

1. The Basic Idea

1.1. The Need for a Predictor

Alice: We now have a code for constant time steps, and one for variable but shared time steps. It's time to bite the bullet, and start working on a code where each particle will have its own individual time step.

This will be quite complicated, but it is certainly necessary. When two particles approach each other closely, we don't want to force all other particles that are far away to take the same small time steps as the closely approaching particles will take.

I've never coded up such a scheme, and I've never seen it in text books on numerical integration either. It seems to be a special trick that is used in astrophysics, in the stellar dynamics community. I presume that you have some experience with . . . why are you smiling?

Bob: Yesterday, after we made such strides with our shared time step code, I just couldn't stop. I've been hacking away the whole evening, and I must admit, a good deal of the night. But I now have a working individual time step version! It's right here under our fingertips, in file nbody_ind1.rb.

Alice: So that's why you are smiling, you have every right to do so!

Bob: And I must admit, yes, I do have some experience. In no way could I have come up with all the tricks of the trade in just one night. The most tricky thing is to visualize the scheduling requirements. But why talk about it abstractely if I can show you concretely what I did?

Alice: Please do!

Bob: First of all, I limited myself to the Hermite scheme.

Alice: Why? After you had derived all those nice integrators, why not carry them over?

Bob: I can see that you don't have experience with individual time step schemes. And in that light, it certainly is a fair question. The answer is that not every algorithm comes with a predictor.

1.2. Crossing the Road

Alice: A predictor?

Bob: Yes. A predictor. Okay, lets start at the basics. If particles have individual time steps, how do you think a single particle will make its next step?

Alice: That will depend on the algorithm.

Bob: Let's say that you have a Runge-Kutta scheme. For definiteness, let us take the second-order Runge-Kutta scheme rk2 that we have been using.

Alice: Well, in that case it will first make a half-step, with a length of .

Bob: How will our particle take that step?

Alice: Is this a riddle?

Bob: No, I'm serious. How does the particle step forward?

Alice: How does the chicken cross the road. Well, of course, you take the velocity and acceleration . . .

Bob: . . . How do you obtain the acceleration?

Alice: In the usual way, as the inverse square of the distance between . . .

Bob: . . . How do you determine the distance?

Alice: You take the position of the other particle and . . .

Bob: . . . How do you know the position of the other particle?

Alice: You look it up! Will you never let me finish a sentence?

Bob: I just did. Now, how do you look it up?

Alice: You look at . . . ahaha! Now I get it. The other particle has a different time step size, so it will have a position all right, but a position that is valid for a different time than the time for which the position of our particle is valid.

Bob: So we have to extrapolate the position of the other particle, to predict where it will be by the time we need its position.

Alice: Hence the need for a predictor.

Bob: Exactly.

Alice: That is tricky indeed. How do you use a Runge-Kutta scheme to provide a predictor?

Bob: The answer is: you don't. Or at least astrophysicists don't. They tend to stick to a very small set of algorithms with known predictors.

Alice: But wait a minute. For a second order Runge-Kutta, it should not be too hard to construct a predictor!

Bob: That may be, but I'm sure it's non-trivial already for a fourth-order Runge-Kutta. And this is why you've never heard of a code for collisional stellar dynamics that uses a Runge-Kutta.

Alice: Ah, is that the reason that people only use Hermite and multi-step methods?

Bob: It's not the only reason, I think, but certainly one of the reasons.

Alice: Not a very good reason, for sure. Just the fact that something is not so easy to do doesn't mean you shouldn't try! I'm happy to start with your Hermite implementation, but I sure would like to try other ones as well. Let's set ourselves a goal, to make a fourth-order Runge-Kutta integrator!

Bob: That's quite a challenge. In the almost half century that astronomers have worked with individual time steps, I've never seen or heard about a Runge-Kutta implementation. But why not? If we succeed, it would mean that in at least one respect we would be ahead of the pack.

Alice: I'd love to see that! But okay, for now let's be content with Hermite. Can you show me what you did?

1.3. A Latest Sales Date

Bob: Let me skip the start-up phase for now. Often in life, it is not so easy to start at the beginning. Let's start in the middle instead. Imagine that all of our particles have already taken at least one step. Let me show you how which particle takes the next step.

Alice: How which what?

Bob: There are two questions here: at any given time we first have to determine which particle needs to take a step, and then we have to figure out how that step can be taken.

Alice: I'm all ears.

Bob: First a bit more background. In the past, the state of each particle was characterized by a mass @mass, position @pos, and velocity @vel, all given at a shared system time @time, an instance variable of the class NBody. In the case of individual time steps, each particle has its own time @time, which is now an instance variable of the class Body.

In addition, since each body has its own time step, it has a predicted time @next_time, again different for each body, which is the time of completion of the next time step: @next_time - time is by definition the indivdual time step of the particle.

The concept of a @next_time is an important one: it plays the role of a latest sales date. If you go to the store to buy some cookies, you can look at the box, and read the latest sales date. If that date is in the past, you probably don't want to buy the cookies, since you can't trust them to be fresh enough.

Alice: Though they won't go bad overnight, as soon as the latest sales date has passed.

Bob: Neither will our particles lose predictability completely, after @next_time, but still, it is a good measure of the time until which we can predict the position of a particle. If another particle needs to know the position of our particle at time t,, there is no problem as long as t <= @next_time.

To sum up: each particle carries with it a predictor that allows it to predict its own position during the interval [ @time, @next_time ].

Alice: If another particle asks for the position of our particle at time t, within that interval, how does our particle communicate its predicted position?

Bob: Through its instance variable @pred_pos. Here is the hand shaking mechanism. First the other particle sends the time for which it wants to get the information. If the time is outside the proper interval, our particle will complain. If the time is within the interval, our particle will oblige, and predict its position and velocity, and store those values in the variables @pred_pos and @pred_vel. The other particle can then read off the values, and determine the acceleration and jerk exercised on itself by our particle.

Alice: Ah yes, jerk, we are dealing with a Hermite scheme, and therefore you have to predict the velocity as well.

Bob: Indeed. For any other scheme that we have implemented so far, it would suffice just to predict the position, since we would only need the acceleration.

Alice: So our Body class now has four extra instance variables, over and above what we had in the shared time stap class: time, @next_time, @pred_pos, and @pred_vel.

Bob: Indeed.

1.4. Scheduling

Alice: You mentioned that you wanted to explain how to find the particle that needs to be moved, before explaining how to move it.

Bob: Yes, but we need still more background first. Let us again call t the time at which we decide to look at the system. We know for sure that for all particles, @time <= t and @next_time >= t.

Alice: Why?

Bob: Whenever a particle runs out of predictive power, it has to be updated immediately. In other words, as soon as the latest sales date, @next_time, is passed, the particle position needs to be updated. This means that the particle takes a step at this time, and in doing so, it acquires a new latest sales date. So if all marches well, the latest sales date for each particle will remain in the future, or at worst will be in the present, never in the past. This means @next_time >= t.

Alice: Fair enough.

Bob: And since you cannot predict particles past their latest sales date, and since latest sales dates can be arbitrarily close to the present time t, no particle can take a step in the future. Only when the time t has caught up with the latest sales date @next_time of a particular particle, will that particle make a step. When it does so, its new time time will take on the same value that @next_time had, which at the time of the step taking is exactly t. From that time on, t will march forward, while @time is frozen until the next step. Hence @time <= t. QED.

Alice: Hmmm. I sort-of get it, but I must admit, I don't see it very clearly yet. I bet looking at the code will makes things clearer.

Bob: It always does. Here is how it works. Let's jump right into the heart of the matter, the method evolve. For comparison, here is what evolve looked like in the case of shared time steps: 

   def evolve(c)
     @nsteps = 0
     @e0 = ekin + epot
     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_flag
 
     while @time < t_end
       @dt = c.dt_param * collision_time_scale
       if c.exact_time_flag and @time + @dt > t_out
         @dt = t_out - @time
       end
       send(c.method)
       @time += @dt
       @nsteps += 1
       if @time >= t_dia
         write_diagnostics
         t_dia += c.dt_dia
       end
       if @time >= t_out - 1.0/VERY_LARGE_NUMBER
         acs_write
         t_out += c.dt_out
       end
     end
   end

And here is how I rewrote it for the individual time step case: 

   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

1.5. Back to How and Which

Alice: What does startup do?

Bob: Let's look at that later. It records the initial total energy @e0@, as in the shared time step case, and it does a lot more. Other than that, the parts up to the while loop are the same in both case. So let us jump to the while loop, for now.

In the case of shared time steps, we first determined the size of the next (shared) time step, and then we pushed all particles forward, through the command send(c.method). Following that, there are only some administrative issues, related to reporting diagnostics and orchestrating output. These issues are quite similar in both cases. Note that in the case of individual time steps, I've made sure to synchronize the particles will a call to sync; but again, let's leave that for later.

Alice: Are you leaving anything for now?

Bob: Yes. Do you remember the how which question?

Alice: You were going to tell me how to push which particle forward. That was quite a while ago.

Bob: Well, we needed to lay some foundations. Now I can show you. The which part is decided right at the start of the while loop, by a call to find_next_particle. The how part is decided three lines further on: the particle np that is found is being asked to take a time step, through its instance method autonomous_step.

Alice: Why autonomous?

Bob: This is the generic situation, well after the start and well before the finish of a calculation. As you can see from the if condition one line earlier, an autonomous step is only allowed if the end of the step np.next_time is still earlier than the finishing time t_end. If that is not the case, we have to wrap it, and start synchronizing everything.

Alice: Can you show me what happens in the generic case, with an autonomous step?

Bob: For that, we have to go to the Body class.
Previous ToC Up Next