Previous ToC Up Next

4. The Answer

4.1. Inspection

Alice: Here is a cup of tea. I figured you'd need one!

Bob: Thanks! I've been scratching my head, but pure thought hasn't given me an answer yet.

Alice: Well, debugging is part of life, at least if your life is tied up with simulations. When in doubt, sort it out. Let's have a look at what is going on under the hood. How about just printing out some output files?

Bob: When in doubt, print it out, you mean. Okay, here are the first two leapfrog files:

 |gravity> cat leap_sh.out
 ACS
   NBody 
     Array body
       Body body[0]
         Fixnum body_id
           1
         Float mass
              2.5000000000000000e-01
         Vector pos
              5.5614242998122632e-01   6.1008593316506865e-01   3.6262238208463954e-01
         Vector vel
             -2.9224887123472459e-01   7.2622108865256441e-01  -8.7148985224906009e-02
       Body body[1]
         Fixnum body_id
           2
         Float mass
              2.5000000000000000e-01
         Vector pos
             -7.2096410104539543e-02   4.3661606633070826e-01   1.1923039259402968e-01
         Vector vel
              6.6764942540701733e-01  -7.2481737546260958e-02   8.7603134996014700e-02
       Body body[2]
         Fixnum body_id
           3
         Float mass
              2.5000000000000000e-01
         Vector pos
              1.2449593150356981e-01   4.9538247509110966e-01   1.8426101006668336e-01
         Vector vel
             -1.2801484357979337e-01  -9.7767564517408012e-01  -4.9766334605570267e-01
       Body body[3]
         Fixnum body_id
           4
         Float mass
              2.5000000000000000e-01
         Vector pos
             -6.0854195138025691e-01  -1.5420844745868865e+00  -6.6611378474535232e-01
         Vector vel
             -2.4738571059249900e-01   3.2393629406777680e-01   4.9720919628459398e-01
     Float dt
          1.4283994962464802e-03
     Float e0
         -2.4999999999999989e-01
     Fixnum nsteps
       54
     Float time
          1.0135606932368554e-01
 SCA
 |gravity> cat leap_sh_ten.out
 ACS
   NBody 
     Array body
       Body body[0]
         Fixnum body_id
           1
         Float mass
              2.5000000000000000e-01
         Vector pos
              5.5652480881381716e-01   6.0913185798706937e-01   3.6273619751754760e-01
         Vector vel
             -2.9034161655073354e-01   7.2676302364413203e-01  -8.6367972342313720e-02
       Body body[1]
         Fixnum body_id
           2
         Float mass
              2.5000000000000000e-01
         Vector pos
             -7.2967298036620262e-02   4.3671275945323129e-01   1.1911719477854188e-01
         Vector vel
              6.6062456540744963e-01  -7.4520519429431789e-02   8.5258939613185855e-02
       Body body[2]
         Fixnum body_id
           3
         Float mass
              2.5000000000000000e-01
         Vector pos
              1.2465960694350207e-01   4.9666505124726201e-01   1.8491313863703351e-01
         Vector vel
             -1.2284109356974235e-01  -9.7602658875333781e-01  -4.9603558166343786e-01
       Body body[3]
         Fixnum body_id
           4
         Float mass
              2.5000000000000000e-01
         Vector pos
             -6.0821711772069864e-01  -1.5425096686875648e+00  -6.6676653093312255e-01
         Vector vel
             -2.4744185528697535e-01   3.2378408453863577e-01   4.9714461439256524e-01
     Float dt
          1.4286878949393660e-04
     Float e0
         -2.4999999999999989e-01
     Fixnum nsteps
       533
     Float time
          1.0004316391104366e-01
 SCA

4.2. Of Course!

Alice: And look, there is the answer! The final output times are quite different.

Bob: Indeed! And of course!! Remember our discussion when we wrote the shared timestep code? We realized that we would overshoot, but we decided not to care. Well, a constant time step scheme can stop properly after one time unit, but a shared time step code is guaranteed to overshoot, unless we teach it to halt at exactly the required time.

Alice: And we decided not to do that, because we did not want to influence the integration. If you restart a code from a previous output, say at time t=1, it would be nice to get the same results at a later time, t=2 for example, as if you had made a single run directly to t=2.

Bob: Yes, that was a good argument, but we can't have our cake and eat it. At least for testing purposes, we should be able to ask the code to stop at exactly the time we want it to stop.

Alice: Perhaps we can have our cake and eat it, if we build in an extra option. Running the code with the option --exact_time, if we want to call it that way, should produce all diagnostics at the exact time that we order them. For our current purposes, that would be ideal. The default should be to have only approximate obedience to time requests, in order to allow accurate restarts. So this extra option should be a boolean flag.

Bob: That's easy to implement. How about this:

   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

I have added an extra if statement at the top of the while loop, which causes the time step to shorten whenever it threatens to overshoot the next output time t_out. Note that I subtracted a very small amount from t_out in the last if statement, just in case round-off errors would make @time ever so slightly smaller than the target time t_out. In this way, even with a small roundoff error, @time should still be larger than t_out - 1.0/VERY_LARGE_NUMBER.

And to make it complete, we should implement the option. Here is the full new list:

 |gravity> kali nbody_sh1.rb -h
 Shared Time Step Code
 -g  --integration_method: Integration method            [default: hermite]
 -c  --step_size_control        : Parameter to determine time step size [default: 0.01]
 -d  --diagnostics_interval: Interval between diagnostics output [default: 1]
 -o  --output_interval  : Time interval between snapshot output [default: 1]
 -t  --duration         : Duration of the integration    [default: 10]
 --exact_time           : Force all outputs to occur at the exact times
 -i  --init_out         : Output the initial snapshot
 --verbosity            : Screen Output Verbosity Level  [default: 1]
 --acs_verbosity                : ACS Output Verbosity Level     [default: 1]
 --precision            : Floating point precision       [default: 16]
 --indentation          : Incremental indentation        [default: 2]
 -h  --help             : Help facility
 ---help                        : Program description (the header part of --help)

4.3. The Right Thing

Alice: Let's inspect an output file, to make sure it halts correctly.

Bob: Let's go back to our original Runge-Kutta tests. Here you are:

 |gravity> kali nbody_sh1.rb -g rk4 --exact_time -t 0.1 -d 0.1 -o 0.1 < test.in > rk4sh.out
 ==> Shared Time Step Code <==
 Integration method             : method = rk4
 Parameter to determine time step size: dt_param = 0.01
 Interval between diagnostics output: dt_dia = 0.1
 Time interval between snapshot output: dt_out = 0.1
 Duration of the integration    : t = 0.1
 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 = 0.1, after 54 steps :
   E_kin = 0.337 , E_pot =  -0.587 , E_tot = -0.25
              E_tot - E_init = 1.05e-10
   (E_tot - E_init) / E_init = -4.19e-10
 |gravity> cat rk4sh.out
 ACS
   NBody 
     Array body
       Body body[0]
         Vector a0
             -1.4551113677428160e+00  -4.0937248040602131e-01  -5.9484396408668083e-01
         Vector a1
             -1.4550361621604517e+00  -4.0952537438487502e-01  -5.9485662726960709e-01
         Vector a2
             -1.4549609178192939e+00  -4.0967823944669401e-01  -5.9486926744033863e-01
         Fixnum body_id
           1
         Float mass
              2.5000000000000000e-01
         Vector old_pos
              5.5655846997385572e-01   6.0904756824642325e-01   3.6274620891188408e-01
         Vector pos
              5.5653733908597036e-01   6.0910048509160519e-01   3.6273992404177713e-01
         Vector vel
             -2.9027882725502885e-01   7.2678071310119674e-01  -8.6342298172199225e-02
       Body body[1]
         Vector a0
              5.2980641547375358e+00   1.5524091751430418e+00   1.7755933059108742e+00
         Vector a1
              5.3001727066872686e+00   1.5524850803288315e+00   1.7760169226324680e+00
         Vector a2
              5.3022830029452255e+00   1.5525608408059413e+00   1.7764407959724633e+00
         Fixnum body_id
           2
         Float mass
              2.5000000000000000e-01
         Vector old_pos
             -7.3043866243215211e-02   4.3672141412104487e-01   1.1910732175994151e-01
         Vector pos
             -7.2995797933024742e-02   4.3671597940599116e-01   1.1911351903693274e-01
         Vector vel
              6.6039554961852320e-01  -7.4587522216318800e-02   8.5182239474103624e-02
       Body body[2]
         Vector a0
             -3.8856720863240786e+00  -1.2589212191004651e+00  -1.2299390761410594e+00
         Vector a1
             -3.8878578159544785e+00  -1.2588464839465736e+00  -1.2303500498675188e+00
         Vector a2
             -3.8900453288247379e+00  -1.2587716330493681e+00  -1.2307613032151099e+00
         Fixnum body_id
           3
         Float mass
              2.5000000000000000e-01
         Vector old_pos
              1.2467381726519004e-01   4.9677823550101824e-01   1.8497065495254192e-01
         Vector pos
              1.2466489594234119e-01   4.9670717983278484e-01   1.8493454651388136e-01
         Vector vel
             -1.2267302289092370e-01  -9.7597227312898527e-01  -4.9598243247223878e-01
       Body body[3]
         Vector a0
              4.2719299329358869e-02   1.1588452436344483e-01   4.9189734316866303e-02
         Vector a1
              4.2721271427661349e-02   1.1588677800261732e-01   4.9189754504657977e-02
         Vector a2
              4.2723243698806396e-02   1.1588903169012071e-01   4.9189774682985216e-02
         Fixnum body_id
           4
         Float mass
              2.5000000000000000e-01
         Vector old_pos
             -6.0818842099583037e-01  -1.5425472178684871e+00  -6.6682418562436785e-01
         Vector pos
             -6.0820643709528655e-01  -1.5425236443303820e+00  -6.6678798959259156e-01
         Vector vel
             -2.4744369947257061e-01   3.2377908224410784e-01   4.9714249117033416e-01
     Float dt
          7.2808425749776307e-05
     Float e0
         -2.4999999999999989e-01
     Fixnum nsteps
       54
     Float time
          1.0000000000000001e-01
 SCA
Alice: Good! Indeed at time 1. Let's repeat our comparison with the constant time step run.

Bob: That should come out a whole lot better now! Here we go:

 |gravity> cat rk4sh.out rk4cst.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.9973298249408804e-10
Alice: And right you are. That's more like it. Wonderful!

4.4. Setting up the Leapfrog

Bob: Let's repeat our leapfrog test suite, just to make sure that everything behaves the way it should. It is hard to use a sixth-order integrator, since before you know it you are up to the round-off barrier.

Alice: I agree. Always good to make an extra test.

Bob: Here are the three output files, now finished at the exact time:

 |gravity> kali nbody_sh1.rb -g leapfrog --exact_time -t 0.1 -d 0.1 -o 0.1 < test.in > leap_sh.out
 ==> Shared Time Step Code <==
 Integration method             : method = leapfrog
 Parameter to determine time step size: dt_param = 0.01
 Interval between diagnostics output: dt_dia = 0.1
 Time interval between snapshot output: dt_out = 0.1
 Duration of the integration    : t = 0.1
 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 = 0.1, after 54 steps :
   E_kin = 0.337 , E_pot =  -0.587 , E_tot = -0.25
              E_tot - E_init = 4.23e-06
   (E_tot - E_init) / E_init = -1.69e-05
 |gravity> kali nbody_sh1.rb -g leapfrog --exact_time -c 0.001 -t 0.1 -d 0.1 -o 0.1 < test.in > leap_sh_ten.out
 ==> Shared Time Step Code <==
 Integration method             : method = leapfrog
 Parameter to determine time step size: dt_param = 0.001
 Interval between diagnostics output: dt_dia = 0.1
 Time interval between snapshot output: dt_out = 0.1
 Duration of the integration    : t = 0.1
 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 = 0.1, after 533 steps :
   E_kin = 0.337 , E_pot =  -0.587 , E_tot = -0.25
              E_tot - E_init = 4.2e-08
   (E_tot - E_init) / E_init = -1.68e-07
 |gravity> kali nbody_sh1.rb -g leapfrog --exact_time -c 0.0001 -t 0.1 -d 0.1 -o 0.1 < test.in > leap_sh_hundred.out
 ==> Shared Time Step Code <==
 Integration method             : method = leapfrog
 Parameter to determine time step size: dt_param = 0.0001
 Interval between diagnostics output: dt_dia = 0.1
 Time interval between snapshot output: dt_out = 0.1
 Duration of the integration    : t = 0.1
 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 = 0.1, after 5330 steps :
   E_kin = 0.337 , E_pot =  -0.587 , E_tot = -0.25
              E_tot - E_init = 4.2e-10
   (E_tot - E_init) / E_init = -1.68e-09

4.5. Testing the Leapfrog

Alice: I bet this will come out correctly quadratic.

Bob: I think so to, but seeing is believing:

 |gravity> cat leap_sh.out leap_sh_hundred.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.4645881582082328e-05
 |gravity> cat leap_sh_ten.out leap_sh_hundred.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.4383516282539215e-07
Alice: We see, we believe, we are delighted!

Bob: You may say so! I am really relieved. I was so frustrated, when we found that linear behavior! But now all is well.
Previous ToC Up Next