Previous ToC Up Next

4. Spacetime

4.1. Mileage

Bob: Yesterday you mentioned that you could see mileage in taking a more abstract, four-dimensional perspective on the N-body problem.

Alice: And you looked skeptical.

Bob: Because I am. What can be the mileage in such a shift, given that we are not dealing with relativistic velocities?

Alice: For one thing, a spacetime view puts our discussion of restarting a run from an output file in a clearer perspective. If we create a new class Worldline for each world line, we can rename our current class Body to Worldpoint and instead of NBody we can talk about World, in order to point to the whole system.

For quite a while, I've been unhappy with the big gap between the bottom-level Body objects and the whole NBody system, as a single top-level objects. With a Worldline class in between, we can separate many issues in a more fine-grained way.

Look at it this way. If you start in four dimensions, identities of objects are related to world lines, while specific values of positions and velocities are related to worldpoints. But as soon as you project that spacetime view into a single space view, you wind up with point particles that have position and velocity and identity, and these three all look the same: the original distinction between lines and points is lost.

In the three dimensional view, we are dealing with moving particles. Each particle has an identity that has to be carried with it while its position and velocity are changing. And with individual time steps, it becomes rather complicated and potentially very confusing, to keep track of which particle is where at what time, and to make sure that particles interact with each other in the right way.

Bob: Can you be more concrete?

Alice: In your individual time step code nbody_ind1.rb, the main engine evolve has to juggle the regular phase of evolving particles with the occasional interruptions of output demand, which require synchronizing all particles. As we have seen, this means that a change in diagnostics frequency has physical repercussions: it means that our particles will be forced to move on slightly different orbits.

This is a Bad Thing.

It would be much better to implement a form of diagnostics reporting that is totally decoupled from the way particles move. If we can make orbit integration and diagnostics reporting truly modular, without the one influencing the other, I would be much happier.

4.2. Eras and Snapshots

Bob: And introducing the concept of a Worldline could do that, you think?

Alice: I think so. We could leave all diagnostics, together with all other output issues to the top level World class. The evolution could take place at a lower level.

Bob: At the Worldline level?

Alice: No, a Worldline class contains only information about a single particle. We need yet again something intermediate, this time in between Worldline and World. Something like a block of time during which all particles can move forward without worrying about output issues. If we save all the information about all steps that are taking during this time, in other words if we save all Worldpoint information, we can then deal with output issues at the end of that block of time.

Bob: Putting a Worldblock class in between World and Worldpoint?

Alice: Yes, but we need a better name. How about a period?

Bob: A dynasty?

Alice: Effectively, yes. It is rather recent that people are counting time in terms of a universal and continuous year count. Most cultures used to talk in terms of the fourth year of the reign of king so-and-so.

Bob: But Worlddynasty doesn't roll of the tongue.

Alice: How about calling it an era? That's about the shortest term I can think of.

Bob: Can't make it much shorter. A Worldera class?

Alice: That sounds better already. So instead of dealing with two levels, of Body and NBody we will now have four levels: Worldpoint, Worldline, Worldera, and finally just World.

Bob: But what about output? I don't see any natural place yet to synchronize particles in this hierarchy.

Alice: Good point. The four world-related classes are all intrinsically temporal constructs. We need a separate class for spatial information. This class should contain the results of taking a hypersurface cut through a bunch of world lines.

Bob: You mean, the result of taking a snapshot of the system?

Alice: Yes, that's a good name! Let's call it a Worldsnapshot class. By using the word `world' everywhere, we make it clear that we're dealing with a coherent system of concepts.

Bob: Before you introduce even more classes, shall we start with a toy code, to see whether your grand vision makes any sense, in practice?

Alice: I'd love to! It shouldn't be too hard to rewrite your nbody_ind1.rb into World form.

Bob: I'm not at all convinced that going from two to five classes will simplify life, but I must admit, your ideas sound intriguing enough to try it out.

Before getting started, let's make sure this time that we obey Ruby's convention of using MixedCaseNotation from the beginning: let's call our classes World, WorldSnapshot, WorldEra, WorldLine, and WorldPoint.

Alice: MyPleasure.

4.3. Intermezzo

From this point on, Alice and Bob spent quite a bit of time coding up an individual time step algorithm, using their five new classes. As one can imagine, they ran into all kinds of problems, and went through many a debug session. Finally, after several days, it all worked. We will join them now, while they look at their first finished product, in file world1.rb. This will be the first of a series of such implementations, and the detailed layout will change quite a bit. However, the current code contains all the basic ingredients for a worldline-based N-body code.

Bob: That was quite a lot of work. The new file world1.rb is a complete overhaul of what we had in nbody_ind1.rb. It is twice a long, too.

Alice: But it is a lot more structured, and I'm sure it will be much easier to extend to more complicated situations. Besides, it has significantly more functionality.

Bob: Before we congratulate ourselves too much, we should make sure that we get the same answers as before. Let us apply our standard test again, using the by now familiar initial conditions for our 4-particle system:

 |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
 ==> 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
              actual seed used  : 1
Here is our comparison run for the older code:

 |gravity> kali nbody_ind1.rb -t 1 -c 0.001 < test.in > tmpi.out
 ==> Individual Time Step Hermite Code <==
 Parameter to determine time step size: dt_param = 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 42565 steps :
   E_kin = 0.313 , E_pot =  -0.563 , E_tot = -0.25
              E_tot - E_init = -1.42e-12
   (E_tot - E_init) / E_init = 5.68e-12
And here is what our new code gives:

 |gravity> kali world1.rb -t 1 -c 0.001 < test.in > tmpw.out
 ==> Individual Time Step, Individual Integration Scheme Code <==
 Determines the time step size  : dt_param = 0.001
 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 42578 steps to time 1):
     E_kin = 0.313 ,     E_pot =  -0.563 ,      E_tot = -0.25
        E_tot - E_init = -1.05e-12
         (E_tot - E_init) / E_init = 4.21e-12
Aha, at least the number of time steps is very similar, and so is the total energy error. Here is the distance in phase space

 |gravity> cat tmpi.out tmpw.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:   4.2386924923927053e-12

4.4. Almost Perfect Scaling

Alice: I'm glad to see that: the distance is comparable to the energy error. In fact, the correspondance is so good, that I wonder whether we've been a bit lucky. I'd feel more comfortable if we would redo the run with somewhat lower accuracy, just to check whether the phase space distance will become larger.

Bob: Why not. When we take time steps that are ten times longer, we get:

 |gravity> kali nbody_ind1.rb -t 1 -c 0.01 < 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
I bet the new code will show a similar error:

 |gravity> kali world1.rb -t 1 -c 0.01 < test.in > tmpw.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
Alice: And so it does. So now let us see whether the phase space difference is also be comparable:

 |gravity> cat tmpi.out tmpw.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:   1.0927259036948475e-08
Bob: Which it is. Okay, we seem to have done something right.
Previous ToC Up Next