Previous ToC Up Next

5. Reproducibility Lost

5.1. Making Sure

Alice: Now that everything seems to work, I would like to get back to this question of reproducibility. We have argued that we can integrate a system for a while, get an output, and then restart a second integration from that ouput to a finishing time, while getting the exact same result as when we would have done a run from the beginning directly up to this finishing time.

Bob: Unless we set the --exact_time flag, in which case this no longer holds.

Alice: I think you're right, but I'd like to make sure. First, I would like to check whether we get exact conservation if we don't set the --exact_time flag. And then I'd like to see how much difference it makes, if we set that flag.

Bob: I'm game. It will certainly give us more confidence in our codes, if they do what they should be doing. But let's start at the beginning, with the constant time step code. Let's run it twice for one time unit, piping the result of the first one into the second one, and once for two time units. We can then compare the results with our phase space measuring rod.

We may as well take the same case we tested before, starting from the file test.in, containing a Plummer model with five particles: 

 |gravity> kali nbody_cst1.rb -t 1 < test.in > test01.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
Now I'll run it for another time unit, till time 2: 

 |gravity> kali nbody_cst1.rb -t 2 < test01.out > test12.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 = 1, after 0 steps :
   E_kin = 0.248 , E_pot =  -0.613 , E_tot = -0.365
              E_tot - E_init = 0
   (E_tot - E_init) / E_init = -0
 at time t = 2, after 1000 steps :
   E_kin = 1.3 , E_pot =  -2.12 , E_tot = -0.815
              E_tot - E_init = -0.45
   (E_tot - E_init) / E_init = 1.23
 at time t = 3, after 2000 steps :
   E_kin = 1.19 , E_pot =  -2.7 , E_tot = -1.52
              E_tot - E_init = -1.15
   (E_tot - E_init) / E_init = 3.15
Hey, that is strange, it ran all the way to time 3, even though I gave it the option -t 2.

5.2. The Danger of Short Options

Alice: That shows the danger of using short options. I bet we defined the argument for this option to be the time difference, not the target time of the integration. Let's check: 

 |gravity> kali nbody_cst1.rb --help -t
 -t  --duration         : Duration of the integration    [default: 10]
 
     This option sets the duration t of the integration, the time period
     after which the integration will halt.  If the initial snapshot is
     marked to be at time t_init, the integration will halt at time
     t_final = t_init + t.
Bob: I really like our help facility. Helpful for sure! Okay, I'll try again: 

 |gravity> kali nbody_cst1.rb -t 1 < test01.out > test12.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 = 1, after 0 steps :
   E_kin = 0.248 , E_pot =  -0.613 , E_tot = -0.365
              E_tot - E_init = 0
   (E_tot - E_init) / E_init = -0
 at time t = 2, after 1000 steps :
   E_kin = 1.3 , E_pot =  -2.12 , E_tot = -0.815
              E_tot - E_init = -0.45
   (E_tot - E_init) / E_init = 1.23
And now I'll run the whole show right through from the beginning: 

 |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
Now if everything behaves as I expect it will, we should get zero distance in phase space: 

 |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

5.3. Know Your Output File

Alice: Not quite.

Bob: That's an understatement! The difference is huge. What happened?

Alice: No idea. But as we did before, we may have to inspect the files again. First, let's see what we have here: 

 |gravity> ll test??.out
 ll: コマンドが見つかりません.
Bob: Ah, the file test02.out is twice as big as the others. Of course! By default, the integrator gives one full output per time unit. This means that the first file, the one we are comparing with by default, will be the result of an integration to time 1, not 2. That explains!

Alice: If it does, we should at least get a null result in a comparison with the first shorter run: 

 |gravity> cat test01.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
Bob: And yes, that's obviously the correct explanation. Well, let me repeat the longer run, but this time with only one output at the end. 

 |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
Wanna bet we'll get zero now?

Alice: I won't hold my breath. Let's check: 

 |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
Bob: You could have held your breath, and you wouldn't have suffered! Good, so we have reproducibility, at least for the constant time step scheme.

5.4. Synchronization

Alice: Let's do a similar thing for shared time steps, first with our exact_time option. 

 |gravity> kali nbody_sh1.rb -t 1 --exact_time < test.in > test01e.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 -t 1 --exact_time < test01e.out > test12e.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 = 1, after 0 steps :
   E_kin = 0.313 , E_pot =  -0.563 , E_tot = -0.25
              E_tot - E_init = 0
   (E_tot - E_init) / E_init = -0
 at time t = 2, after 1224 steps :
   E_kin = 0.171 , E_pot =  -0.421 , E_tot = -0.25
              E_tot - E_init = 2.86e-10
   (E_tot - E_init) / E_init = -1.14e-09
 |gravity> kali nbody_sh1.rb -t 2 -o 2 --exact_time < test.in > test02e.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
 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.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, after 3151 steps :
   E_kin = 0.171 , E_pot =  -0.421 , E_tot = -0.25
              E_tot - E_init = 5.78e-11
   (E_tot - E_init) / E_init = -2.31e-10

5.5. The Toll it Takes

Bob: According to our predictions, the composition of the first two shorter runs should not give exactly the same result as the longer runs. Let's check: 

 |gravity> cat test12e.out test02e.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:   1.3922126211172806e-10
Indeed, the difference is slight, but clearly noticeable.

Alice: Let's make sure that it's not a matter of having a different ending time: 

 |gravity> cat test12e.out test02e.out
 ACS
   NBody 
     Array body
       Body body[0]
         Fixnum body_id
           1
         Float mass
              2.5000000000000000e-01
         Vector old_acc
              3.6193191054803503e-01  -1.9760997053431525e-01   1.9221511129758645e-01
         Vector old_jerk
             -1.6132599715142393e-01   3.4113042593532333e-01   2.6538251242837169e-01
         Vector old_pos
             -3.7432412047693525e-01   3.2999341450049463e-01  -3.0728986598861008e-01
         Vector old_vel
             -1.7759637111914583e-01  -6.5838867169294302e-01  -3.3323751222978998e-01
         Vector pos
             -3.7473997754816485e-01   3.2844750276597751e-01  -3.0807151012094697e-01
         Vector vel
             -1.7674729095980554e-01  -6.5885155999717571e-01  -3.3278561416074259e-01
       Body body[1]
         Fixnum body_id
           2
         Float mass
              2.5000000000000000e-01
         Vector old_acc
              1.3771842548785003e+00   1.1182787946202655e+00   7.3859480075224715e-01
         Vector old_jerk
              4.0921523631915013e+00  -4.3681144401613081e-02   3.4327676374029420e-01
         Vector old_pos
              5.0980674567049045e-01   8.1814032270369028e-02  -7.4105774591426893e-02
         Vector old_vel
              1.0210186743120495e-01   2.0229197405231467e-01   2.4598136277311954e-02
         Vector pos
              5.1005020147075963e-01   8.2291931865753523e-02  -7.4046002568255234e-02
         Vector vel
              1.0534577743505329e-01   2.0491671781238224e-01   2.6332745064625546e-02
       Body body[2]
         Fixnum body_id
           3
         Float mass
              2.5000000000000000e-01
         Vector old_acc
             -1.9365553672199651e+00  -1.1062823001550299e+00  -8.2544760102850889e-01
         Vector old_jerk
             -4.0220898817748783e+00  -2.6958701878300001e-01  -3.5372867640179684e-01
         Vector old_pos
              7.7686658629667216e-01   2.5093926149359153e-01   4.8066063082454383e-02
         Vector old_vel
              9.9818279483534625e-02  -1.6786776488798624e-01  -2.0490858593868053e-01
         Vector pos
              7.7709553619284533e-01   2.5054219468031019e-01   4.7582827583393508e-02
         Vector vel
              9.5261609328486291e-02  -1.7046521400616452e-01  -2.0684708396688911e-01
       Body body[3]
         Fixnum body_id
           4
         Float mass
              2.5000000000000000e-01
         Vector old_acc
              1.9743920179342983e-01   1.8561347606907969e-01  -1.0536231102132472e-01
         Vector old_jerk
              9.1263515734800932e-02  -2.7862262750710284e-02  -2.5493059976686899e-01
         Vector old_pos
             -9.1234920352799009e-01  -6.6274670263674507e-01   3.3332958158275522e-01
         Vector old_vel
             -2.4323767959298123e-02   6.2396446810478179e-01   5.1354796590408236e-01
         Vector pos
             -9.1240575213480857e-01  -6.6128162367124199e-01   3.3453468920040097e-01
         Vector vel
             -2.3860087966902522e-02   6.2440006176754859e-01   5.1329995707622100e-01
     Float dt
          2.3471973126603096e-03
     Float e0
         -2.5000000023295360e-01
     Fixnum nsteps
       1224
     Float time
          2.0000000000000000e+00
 SCA
 ACS
   NBody 
     Array body
       Body body[0]
         Fixnum body_id
           1
         Float mass
              2.5000000000000000e-01
         Vector old_acc
              3.6181397424436307e-01  -1.9736066675347447e-01   1.9240893406399123e-01
         Vector old_jerk
             -1.6168077177383272e-01   3.4166800735759456e-01   2.6546301088387902e-01
         Vector old_pos
             -3.7445371217931045e-01   3.2951257928411237e-01  -3.0753315849144236e-01
         Vector old_vel
             -1.7733211656352585e-01  -6.5853288363726170e-01  -3.3309707806634120e-01
         Vector pos
             -3.7473997754875182e-01   3.2844750276637186e-01  -3.0807151012102341e-01
         Vector vel
             -1.7674729096114672e-01  -6.5885155999636025e-01  -3.3278561416097957e-01
       Body body[1]
         Fixnum body_id
           2
         Float mass
              2.5000000000000000e-01
         Vector old_acc
              1.3801866065485837e+00   1.1182522252318787e+00   7.3884985436282846e-01
         Vector old_jerk
              4.1307974764229956e+00  -2.9083433859158590e-02   3.5527518249497908e-01
         Vector old_pos
              5.0988167212043201e-01   8.1962052359618828e-02  -7.4087615073068186e-02
         Vector old_vel
              1.0310863858202257e-01   2.0310857692795328e-01   2.5137581205126873e-02
         Vector pos
              5.1005020146824909e-01   8.2291931855993844e-02  -7.4046002574291558e-02
         Vector vel
              1.0534577736436623e-01   2.0491671775990278e-01   2.6332745027608583e-02
       Body body[2]
         Fixnum body_id
           3
         Float mass
              2.5000000000000000e-01
         Vector old_acc
             -1.9395064057367308e+00  -1.1064845601023010e+00  -8.2571026134815972e-01
         Vector old_jerk
             -4.0603220119473162e+00  -2.8437083937674112e-01  -3.6565645284204851e-01
         Vector old_pos
              7.7693896112682814e-01   2.5081638256106198e-01   4.7916210283961917e-02
         Vector old_vel
              9.8403051206855149e-02  -1.6867569096918125e-01  -2.0551145714043556e-01
         Vector pos
              7.7709553619817673e-01   2.5054219469183614e-01   4.7582827591061763e-02
         Vector vel
              9.5261609402715164e-02  -1.7046521395238845e-01  -2.0684708392817655e-01
       Body body[3]
         Fixnum body_id
           4
         Float mass
              2.5000000000000000e-01
         Vector old_acc
              1.9750582494378402e-01   1.8559300162389653e-01  -1.0554852707866003e-01
         Vector old_jerk
              9.1205307298153626e-02  -2.8213734121694801e-02  -2.5508174053680960e-01
         Vector old_pos
             -9.1236691309775242e-01  -6.6229100857082768e-01   3.3370456737020027e-01
         Vector old_vel
             -2.4179565386472014e-02   6.2410000325711046e-01   5.1347095801630671e-01
         Vector pos
             -9.1240575213480146e-01  -6.6128162367121590e-01   3.3453468920039642e-01
         Vector vel
             -2.3860087966971380e-02   6.2440006176753271e-01   5.1329995707624976e-01
     Float dt
          1.6169560864454091e-03
     Float e0
         -2.4999999999999989e-01
     Fixnum nsteps
       3151
     Float time
          2.0000000000000000e+00
 SCA
Bob: No, they both end nicely at time 2. So the difference is really a difference between the orbits. This is the toll that synchronisation takes!
Previous ToC Up Next