Previous ToC Up Next

6. Reproducibility Gained

6.1. Asynchronicity

Alice: We have seen that constant time steps allow you to stop and restart, and to get exactly the same results as you would have gotten, had you done one single integration from being to end. But when we tried to do the same thing with adaptive, shared time steps, we found that we could not reproduce the same orbits.

Bob: And we knew why. In fact, we had predicted that we couldn't.

Alice: Yes. Now the question is: can we regain our reproducibility?

Bob: How?

Alice: What did not work was to force an output at exactly the right time. That had the side effect of forcing the last step to be smaller than it would have otherwise been.

However, there is an alternative. How about letting the system overshoot, so that it can do an output a time it is happy with. When we then restart the system, the integration should happen just like it would have otherwise.

Bob: I see what you mean. Well, let's test it!

Alice: Here is what we did before, but this time without the --exact_time option: 

 |gravity> kali nbody_sh1.rb -t 1 < test.in > test01n.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
 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.00035, after 1927 steps :
   E_kin = 0.312 , E_pot =  -0.562 , E_tot = -0.25
              E_tot - E_init = -2.26e-10
   (E_tot - E_init) / E_init = 9.03e-10
 |gravity> kali nbody_sh1.rb -t 1 < test01n.out > test12n.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
 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 = 1.00035, after 0 steps :
   E_kin = 0.312 , E_pot =  -0.562 , E_tot = -0.25
              E_tot - E_init = 0
   (E_tot - E_init) / E_init = -0
 at time t = 2.00117, after 1224 steps :
   E_kin = 0.171 , E_pot =  -0.421 , E_tot = -0.25
              E_tot - E_init = 2.83e-10
   (E_tot - E_init) / E_init = -1.13e-09
 |gravity> kali nbody_sh1.rb -t 2 -o 2 < test.in > test02n.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 = 2.0
 Duration of the integration    : t = 2.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.00035, after 1927 steps :
   E_kin = 0.312 , E_pot =  -0.562 , E_tot = -0.25
              E_tot - E_init = -2.26e-10
   (E_tot - E_init) / E_init = 9.03e-10
 at time t = 2.00117, after 3151 steps :
   E_kin = 0.171 , E_pot =  -0.421 , E_tot = -0.25
              E_tot - E_init = 5.76e-11
   (E_tot - E_init) / E_init = -2.31e-10
Bob: So you think we'll get a smaller distance now?

Alice: As close to zero as we can!

Bob: Well, let's hope for the best: 

 |gravity> cat test12n.out test02n.out | kali nbody_diff.rb
 ==> 6N-dimensional phase space distance between two N-body systems <==
 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

6.2. Overshooting

Alice: Hmmm, that's not good at all. What did we do to deserve that!

Bob: If we look at the diagnostics output, we see that we needed the same number of steps to reach time 1, exactly 1263 in each case. Then in the continued run, we took another 669 steps. This gives us a total number of 1932 steps. Aha! That is not the number of steps that the second run took.

Alice: Good catch! The second run stopped too early. And yes, of course, we could have noticed that directly, had we looked at the final time: the second run indeed shows an earlier ending time than the first one.

Bob: And that makes sense, too: we ask each run to go fro one or two time units, or slightly more if they overshoot. Now asking a code to overshoot twice in a row is likely to produce a larger overshoot than doing it only once.

Alice: Again, you must be right. Going from time 0 to 1, with an overshoot, and then adding one more time unit plus a new overshoot is not the same as going from 0 to 2 immediately, with whatever overshoot that gives.

Bob: Well, shall we try to let the long run proceed a little further, to match up with the previous one?

Alice: You can certainly try. So this means adding the difference of the two ending times to the desired duration; or in fact slightly less, since we'll produce an overshoot anyway.

Bob: Okay, how about this: 

 |gravity> kali nbody_sh1.rb -t 2.0037 -o 2.0037 < test.in > test02n.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 = 2.0037
 Duration of the integration    : t = 2.0037
 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.00035, after 1927 steps :
   E_kin = 0.312 , E_pot =  -0.562 , E_tot = -0.25
              E_tot - E_init = -2.26e-10
   (E_tot - E_init) / E_init = 9.03e-10
 at time t = 2.00117, after 3151 steps :
   E_kin = 0.171 , E_pot =  -0.421 , E_tot = -0.25
              E_tot - E_init = 5.76e-11
   (E_tot - E_init) / E_init = -2.31e-10
 |gravity> cat test12n.out test02n.out | kali nbody_diff.rb
 ==> 6N-dimensional phase space distance between two N-body systems <==
 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:   9.2780296991970546e-03

6.3. Living in a Logical World

Alice: The right distance, this time, namely zero. Perfect!

Bob: It literally couldn't have been better.

Alice: Wait, wait: we still have only 1927 steps, and we still stop at too early a time. What is going on? According to the diagnostics, nothing has changed from the previous run.

Bob: Now that is strange. Obviously, the phase space distance is zero. So something has changed. Yet I see your point. This is a real paradox.

Alice: What could possibly be . . . .

Bob: . . . Ahaha! Of course! We changed the ending time and output time, by a slight amount, but we did not change the diagnostics output time!

Alice: Of course, of course! Ah, that is a relief. Good! Mystery solved, hopefully. You'd better check, though!

Bob: Okay, here we go again: 

 |gravity> kali nbody_sh1.rb -t 2.0037 -d 2.0037 -o 2.0037 < test.in > test02n.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 = 2.0037
 Time interval between snapshot output: dt_out = 2.0037
 Duration of the integration    : t = 2.0037
 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 = 2.00395, after 3152 steps :
   E_kin = 0.171 , E_pot =  -0.421 , E_tot = -0.25
              E_tot - E_init = 5.75e-11
   (E_tot - E_init) / E_init = -2.3e-10
 |gravity> cat test12n.out test02n.out | kali nbody_diff.rb
 ==> 6N-dimensional phase space distance between two N-body systems <==
 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:   9.2780296991970546e-03
Alice: All is well.

Bob: How easy to get confused!

Alice: And how nice to be reassured, afterwards, that the world is still logical, after all.

6.4. An ACS tail Version

Bob: You know, there is a better solution for the problem we encountered, when we had to deal with an output file that had more than one snapshot in it.

Alice: You mean in the case in which we ran a code for two time steps, but we had forgotten that we get by default an output at every time step, so that when we looked at the start of the file, we found the output for time 1, rather than for time 2?

Bob: Yes, that was the problem. I'm sure we'll run into this again, often, and what we really need is a tool that picks up the last ACS object from a file, not the first one. Normally, if you open a file and start reading, you would read in the first ACS file first, but instead, we need to do something like what in Unix the command tail does, showing the end of a file.

Alice: The Unix command tail has a default of showing you ten lines, but you can ask it to show more or fewer lines, with the option -n as in: 

 |gravity> tail -n 3 test02n.out
     Float time
          2.0039487382209993e+00
 SCA
which shows the last three lines.

Bob: Okay, let us call the file acstail.rb. How about this? 

 #!/usr/local/bin/ruby -w
 
 require "acs"
 
 options_text= <<-END
 
   Description: Returns the last n chunks of ACS data
   Long description:
     This program accepts ACS data, in the form of a series of chunks, each
     chunk starting with "ACS" on the first line, and ending with "SCA" on
     the last line.  Like the Unix program `tail', it returns the last n chunks.
 
     (c) 2004, Piet Hut, Jun Makino; see ACS at www.artcompsi.org
 
     example:
                kali #{$0} -n 5
 
 
   Short name:           -n
   Long name:            --number
   Value type:           int
   Default value:        1
   Description:          Number of last ACS data chunks returned
   Variable name:        n
   Long description:
     This option allows the user to choose the number of ACS chunks that
     will be returned by this program, after a series of ACS chunks have been
     read in.
 
 
   END
 
 clop = parse_command_line(options_text)
 
 chunk = []
 a = ACS_IO.acs_read
 while a
   chunk.push(a)
   chunk.shift if chunk.length > clop.n
   a = ACS_IO.acs_read
 end
 while chunk[0]
   chunk.shift.acs_write
 end

6.5. Testing tail

Alice: Looks fine to me. How about doing our previous experiment, without adjusting the output frequency, and using acstail.rb instead?

Bob: Good idea! Here you go. Let's be systematic. We know that this works: 

 |gravity> kali nbody_cst1.rb -t 2 -o 2 < test.in > test02.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 = 2.0
 Duration of the integration    : t = 2.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
 at time t = 2, after 2000 steps :
   E_kin = 1.3 , E_pot =  -2.12 , E_tot = -0.815
              E_tot - E_init = -0.565
   (E_tot - E_init) / E_init = 2.26
 |gravity> cat test12.out test02.out | kali nbody_diff.rb
 ==> 6N-dimensional phase space distance between two N-body systems <==
 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
and we also know that this doesn't work: 

 |gravity> kali nbody_cst1.rb -t 2 < test.in > test02.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 = 2.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
 at time t = 2, after 2000 steps :
   E_kin = 1.3 , E_pot =  -2.12 , E_tot = -0.815
              E_tot - E_init = -0.565
   (E_tot - E_init) / E_init = 2.26
 |gravity> cat test12.out test02.out | kali nbody_diff.rb
 ==> 6N-dimensional phase space distance between two N-body systems <==
 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:   2.7562194398214426e+00
And if all is well, our new tool should be able to repair the last problem, as follows: 

 |gravity> kali nbody_cst1.rb -t 2 < test.in | kali acstail.rb > test02.out
 ==> Returns the last n chunks of ACS data <==
 Number of last ACS data chunks returned: n = 1
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
 ==> 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 = 2.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
 at time t = 2, after 2000 steps :
   E_kin = 1.3 , E_pot =  -2.12 , E_tot = -0.815
              E_tot - E_init = -0.565
   (E_tot - E_init) / E_init = 2.26
 |gravity> cat test12.out test02.out | kali nbody_diff.rb
 ==> 6N-dimensional phase space distance between two N-body systems <==
 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
Alice: Well done!
Previous ToC Up Next