Bob: So far we have seen only the outer wrapping level, in the World class, where IO was scheduled and the command was given to the WorldEra class to do a bunch of work. The World class really behaves like a top executive.

Alice: Or like a secretary, letting people in and out.

Bob: Or a concierge or doorman then. We've got a rather socialistic code, if all these functions are combined in one! Well, let's see what the WorldEra class looks like.

Alice: Let's start with evolve. As we've done before, it's often easier to start in the middle than at the beginning.

Bob: In the generic case, shared_flag will be false. We have added the option -a to ask the code to run with shared time steps, to make it easier to compare the results with our previous code nbody_sh1.rb. But the default value of this flag is false, in which case the unless statement evaluates to being true. The next line then tells the shortest extrapolated world line to extend itself. After that, we do bookkeeping, by adding one to the variable that counts the time steps.

Alice: That line indeed reads like English. Let's reconstruct the meaning. We need to find the first particle that needs to be pushed along its orbit, from its current world point to the next world point. This means that we needt to find the particle with the earliest sales date, which is the smallest time @next_time, an instance variable of class WorldPoint. However, at this level we are still two levels removed from that class, and we certainly don't want to rely on the particular way that this latest sales date is implemented.

Bob: At this point, we'd better remind our students what this latest sales date means.

Alice: It is the latest time for which our extrapolation from the last world point is still valid. What can be extrapolated till this time are the positions and velocities, using an extrapolation method which is again hidden on this level, and implemented on the next level down, in class WorldLine. All we can do on this higher level is to ask a particle . . .

Bob: . . . I thought particles were two levels down!

Alice: You're right! Thanks. I have to learn to be more precise. All we can do on this higher level is to ask a worldline to extend itself. The way it extends itself is by asking the particle as defined at the last world point to take a step to the next world point. So indirectly we will give a command to the world point.

Now to be precise, yes, you are right, the world point is two levels down. But you can also say that a particle, with a particular identity, is represented in our code by a world line. A single world line corresponds to the history of a single particle. So in that sense, a world line in 4-D stands for a particles, just like a body stood for a particle in our 3-D way of coding.

Bob: I'm still wondering whether we've not introduced too many legalities and abstractions here, but I'll reserve judgment. So to sum up, we look for the particle with the earliest latest sales date. Hmm, listening to myself, I have to admit that that does sound a bit contradictory.

Alice: But it is correct. Each particle has a latest sales date. And from among the ensemble of all particles, we choose the particle that has the earliest such date. But I agree, `earliest latest' is confusing. So we could say instead that we are looking for a particle that has the earliest cutoff time for its extrapolation.

Bob: And if you look toward the future, the world line that has the earliest cutoff is shortest. Hence we ask for the shortest extrapolatable world line. I see. But why did we call it extrapolated?

Alice: I guess extrapolatable just sounds a bit funny, not altogether palatable. We can later change the wording, if we want. Anyway, that particle, once we've found it, is asked to extend itself, and we'll see on the WorldLine level just how it wants to do that.

Bob: Quite a long discussion, to explain a single line of code! And especially a line of code that already looks like plain English. It goes to show that, once you understand how a code works, it is easier to read the source itself than the explanation.

Alice: Isn't that true for anything?

Bob: I guess. Let's move on to the else clause. We will enter the next part when the shared_flag is true, that is, when we want to use shared time steps. In that case, let's see, aha, we again find our `earliest latest salesdate' particle, but this time we ask it a more complicated question. We ask it for the next_time of the last world point on its world line. Funny how you have to read Ruby backwards in order to translate it into an English sentence.

Alice: Looking at it now, I must say, I don't like the fact that we are specifically addressing the array worlpoint[] in the class WorldLine. While this gives us access to the last point in the array, we wind up with a object of class WorldPoint and then we ask this object for its next_time. This is precisely the sort of level crossing that we wanted to avoid.

Bob: Well, the most important thing is: the code works, as we have already tested. Let's not touch it for now. I'm already happy if I can get an overview of all that we've done so far.

Alice: Fair enough, but I do want to get back to this blemish, at some point. This is not modular programming!

Bob: For now, all I want to do is to understand what we wrote here. In the case of shared time steps, shouldn't we ask all particles for their prefered time step, and then take the shortest time step, so as to make nobody unhappy?

Alice: Well, if all particles would be synchronized, what you just said would be correct, but in general, they will not be synced. Each particle will have been updated to their own time @time, and each particle is thinking about taking a next step, in the future, to a time @time_next.

So in order to make nobody unhappy, we should avoid making a shared system step that would go beyond anybody's @time_next. This implies that our shared step should exactly land on the earliest value @time_next within the ensemble of all particles. And this is precisely what we are doing here: the particle with the shortest extrapolated world line has by definition the smallest @time_next value.

Bob: I see, yes, that must be right. So the only difference with the statement three lines higher is that now we obtain the actual time @time_next, and then we ask all particles to make a step toward this time.

Alice: What we really should do is to ask all world lines to extend themselves to this time. What we actually are doing here is to cross two levels and use the next_time method to set the variable @time_next within each particle by hand. To be cleaned up!

Bob: Okay, okay, later, as we already said. Here at the end of the shared time step case, we update the step counter by one. Shouldn't that be updating by N, for N particles? Ah, no, the update happens within the loop, and we pass through the loop N times. Fine.

Alice: let's get back to the beginning of our while loop. We have not yet looked at the conditional statement there. The code tells us that we keep looping only while the last world point for the shortest interpolated world line has a time that is earlier than the time @end_time at which the era ends.

Bob: And now you're going to tell us that we're crossing levels again.

Alice: Don't worry, I'll shut up.

Bob: For now.

Alice: For now, yes. Shortest interpolated, what again did that mean? Ah, it must mean the time up till the last computed world point. Up to that point, positions and velocities can be interpolated between two world points; beyond that point, we have to use extrapolation.

Bob: So shortest interpolated means earliest interpolatable, yes?

Alice: I guess we decided to avoid that word, but if you want to use interpolatable, then we are asking for the earliest time at which a particle runs out of being interpolatable. Now that sounds like a truly awful sentence. What we are really doing is asking for the earliest time at which a particle runs out of interpolation. Or simply, we ask for the earliest actually computed world point among those points that are at the end of a world line.

Bob: Aha, the earliest latest computed point!

Alice: Yes, I you like. Now when the time corresponding to this world point is at least as large as the end time @end_time of our era, we can stop. In this case we know for sure that all particles have passed the finish line. What we really need is a guarantee that every worldine is fully interpolatable from start to finish of an era, during the whole time span from @begin_time till @end_time.

Bob: I remember that we had a discussion about that. We first thought that it might be enough to let every particle make just one step to or beyond the finish. But then we realized that some particles would have to step quite a bit beyond the finish.

Alice: Yes. As long as there is even one particle left that has not yet stepped out of the duration of our era, we cannot stop. But it is possible that this last particle has a rather long time step, and hence a very late latest sales date. Given our algorithm, this particle will only step forward when the system time reached its latest sales date. If this date is way beyond the end of our era, all other particles may be forced to take many steps each, before this straggler particle is finally asked to step out of the era.

Bob: In itself, this is not such a bad situation, since all the work we do now, in the current era, is something we don't have to do in the next era. It is only at the end of a calculation, that we do a bit too much work. And in addition, each era will need more storage than strictly necessary. So for both reasons we decided to put an upper limit to the length of a time step for any particle. Ah yes, this is the parameter dt_max that is passed as a parameter to evolve.

Alice: In methods World#setup_from_world and World#setup_from_snapshot, the value of dt_max is set as being equal to the duration of the era, multiplied by a factor dt_max_param. By default, this last number is equal to one, but we can set it with the command line argument -m or --max_timestep_param.

Bob: So we can limit the growth of an era in two ways, by setting the length dt_era itself, which is the duration during which all world lines are interpolatable. In other words, this is a slice of spacetime between two fixed-time hypersurfaces. Each world line will stick out at both sides, back into the past and forward into the future. And we are guaranteed that each world line has at least one world point on its world line sticking out into the past, before the slice, and one world point sticking out into the future, beyond the slice.

In general, many world lines will have several points sticking out at one or two sides of the slice. All we can say is that the longest line sticking out at either side of the slice has a length that is at most a factor dt_max_param longer than the thickness of the slice.

Alice: Yes, that's the picture. Individual time step algorithms sure complicate matters!

Bob: And so does spacetime visualization. But I must admit, all this has given me a fresh way of looking at what I have been doing in the past, when I coded up individual time step methods. And I do think I have a more complete overview of the whole picture now. Whether or not it will turn out to be useful, has to be seen, but it is fun to see an old idea in a new light.

Alice: For the time being, I think the use of dt_max is okay, but I must say, I'm a bit concerned about the way we let the user specify it. It may seem natural to couple the value of dt_max directly to the value of the era length dt_era, making them equal by default, but there is a danger.

Bob: What type of danger?

Alice: Ideally, when you change the length of an era, by using the command line option -e or --era_length, you don't expect the orbits of the stars to change. But in fact, they will change, initially by a small amount, but because of the exponential instability of stellar dynamics, after a while the orbits will be completely different. That is a bad thing.

Bob: I see. Doubling the length of an era, without changing the option -m or --max_timestep_param would lead to an effective doubling of dt_max, which could lead to larger timesteps for some particles, and hence different orbits. The way to counter this would be to make that last option smaller, also by the same factor of two. Let me try it:

 |gravity> kali mkplummer.rb -n 4 -s 1 | kali nbody_set_id.rb > tmp.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
 |gravity> kali world1.rb -t 1 < tmp.in > tmp1.out
 ==> Individual Time Step, Individual Integration Scheme Code <==
 Determines the time step size  : dt_param = 0.01
 Duration of an era             : dt_era = 0.01
 Maximum time step (units dt_era): dt_max_param = 1.0
 Diagnostics output interval    : dt_dia = 1.0
 Snapshot output interval       : 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 (from interpolation after 0 steps to time 0):
     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 (from interpolation after 4291 steps to time 1):
     E_kin = 0.313 ,     E_pot =  -0.563 ,      E_tot = -0.25
        E_tot - E_init = -1.98e-10
         (E_tot - E_init) / E_init = 7.92e-10
 |gravity> kali world1.rb -t 1 -e 0.02 < tmp.in > tmp2.out
 ==> Individual Time Step, Individual Integration Scheme Code <==
 Determines the time step size  : dt_param = 0.01
 Duration of an era             : dt_era = 0.02
 Maximum time step (units dt_era): dt_max_param = 1.0
 Diagnostics output interval    : dt_dia = 1.0
 Snapshot output interval       : 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 (from interpolation after 0 steps to time 0):
     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 (from interpolation after 4268 steps to time 1):
     E_kin = 0.313 ,     E_pot =  -0.563 ,      E_tot = -0.25
        E_tot - E_init = 1.6e-10
         (E_tot - E_init) / E_init = -6.4e-10
 |gravity> kali world1.rb -t 1 -e 0.02 -m 0.5 < tmp.in > tmp3.out
 ==> Individual Time Step, Individual Integration Scheme Code <==
 Determines the time step size  : dt_param = 0.01
 Duration of an era             : dt_era = 0.02
 Maximum time step (units dt_era): dt_max_param = 0.5
 Diagnostics output interval    : dt_dia = 1.0
 Snapshot output interval       : 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 (from interpolation after 0 steps to time 0):
     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 (from interpolation after 4291 steps to time 1):
     E_kin = 0.313 ,     E_pot =  -0.563 ,      E_tot = -0.25
        E_tot - E_init = -1.98e-10
         (E_tot - E_init) / E_init = 7.92e-10
 |gravity> cat tmp1.out tmp2.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.9922800657753722e-09
 |gravity> cat tmp1.out tmp3.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:   0.0000000000000000e+00

Bob: Perhaps. But it does work!

Alice: It works for a factor of two alright, but if you introduce a factor of three increase for the era length, and you want to compensate by shrinking the maximum timestep parameter by a factor three, you may not be able to cancel the two factors precisely without round-off error.

Bob: I agree that we could come up with a better version. However, as you said, for now the current choice of command line arguments will do.

Alice: Let's move to the beginning of the WorldEra class, to see what other methods we have here. The method initialize, which is called when we order a WorldEra.new sets up an empty array @worldline, ready to accept a bunch of world lines.

The methods acc_and_jerk and timescale appear here only as a way to pass them on to a snapshot. When the request comes, from a lower level, to provide acceleration and jerk, or a time scale value, then first a snapshot is taken, and then the snapshot is asked to care care of the requests.

The method startup_and_report_energy is almost a tautology: it hands on the request to start things up down to the WorldLine level, and it takes a snapshot in order to ask that snapshot to provide the total energy value. This value is then returned to the calling function in World, one level higher. I must say, I still have to remind myself that in Ruby you don't have to write something like return energy: you can leave out the return part, and the value that is automatically returned is the value that the last method evaluated provides, in this case the method total_energy.

Bob: But once you get used to it, Ruby does provide a nicely compact notation. It allows you to pack a lot of information on one page, and still in such a way that it is pretty clear what is being done.

We now come to the next two methods, shortest_extrapolated_worldline and shortest_interpolated_worldline. They speak for themselves: they do just what you think they should do.

Alice: At least when you understand the concepts they are pointing too; but we've just gone over that. evolve we've analyzed already. The method next_era creates a new instance of the class WorldEra. This method is invoked, whenever we have evolved a complete era, and we are ready for another one. Now where does this method gets called again?

Bob: At the end of evolve. It is part of what is handed back to the World level, together with the value of the counter of the number of steps that have been taken.

Alice: But next_era seems to contruct, as the name implies, a next era, the next slice in time, after the current era. And it populates that era with those parts of the world lines that fit into the next era, using WorldLine#next_worldline. Now I'm puzzled. If evolve returns an almost empty new time slice, with only the beginning parts of each world line populated, haven't we lost most of what we just computed?

Bob: We have to go back to the World level, to decide that. Let me print out World#evolve:

Bob: The last few methods in WorldEra are all related to IO through the use of snapshots. At the bottom we see how we do input, starting from a snapshot, in setup_from_snapshot. The real work is done on the WorldLine level, in setup_from_single_worldpoint, where we add each next body to a brand new world line, created in the line above. The line below shows how this new world line is than added at the end of our array of world lines, @worldline.

Alice: I keep wondering whether we shouldn't have choosen the plural version of the name, to make it @worldlines instead. After all the array @worldline[] contains many world lines.

Bob: But when we talk about the 3rd world line, it is nice to be able to write @worldline[2]. It would make less sense to write @worldlines[2]. This is the same logic we used way back when when we gave the NBody class an array body for its bodies.

Alice: Yeah, I can see the logic, either way. Oh well, we started off this way, so we may as well continue.

Bob: Let us now turn to output, and the way it is mediated by snapshots. We have the function take_snapshot(time) that constructs a full snapshot at a given time. It translates into using the generic method take_snapshot_except, except that nothing is excepted, because we give the nil argument for the exception, if you see what I mean.

Alice: Hard to hear what you mean, exceptionally so, but yes, I see what you may have meant. And the reason for this overly complex redirection is that the function take_snapshot_except is the real workhorse that is used in preparation for a force calculation: it presents a single particle with the positions of all other particles, after which that single particle can compute the distances to all other particles, and based on that, accelerations and jerks.

Bob: And this ends our tour through WorldEra. Time to descend one more level, to WorldLine.

Previous

ToC

Up

Next