| Previous | ToC | Up | Next |
Bob: It was fun to play with so many different versions, but I’m beginning to get a little confused as to which version did what. Maybe we should pick just one version, and use that while we add more features and move toward specific applications.
Alice: I agree. And since we decided not to worry, for now at least, about performance, let us concentrate on clarity of expression. I must say, I like the last one best, in file rknbody6.rb, where everything fits on one line. However, the version where we gave the Body class explicit helper variables, in file rknbody4.rb, was even shorter.
Bob: Yeah, I did not particularly like the idea of giving this poor Body class all possible helper variables for all possible application. Even though we discussed more clever ways to do this, frankly, I don’t care too much what we will choose in the end. I liked your suggestion to send a command string to be evaluated dynamical, thereby generating the proper helper variables, mostly because of its novelty.
Alice: And it goes with the spirit of the times: just-in-time-delivery! But what I like best about this latest method is that it obeys the DRY principle: we are not repeating ourselves.
Bob: Apart from the fact that you repeatedly bring up particular principles.
Alice: I’ll repeatedly ignore that. Apart from that point, when clarity is really the criterion, I am not sure whether the last version is really clearer. Let us compare the forward Euler algorithm in both cases, in rknbody4.rb:
def forward(dt)
@body.each{|b| b.old_acc = b.acc(@body)}
@body.each{|b| b.pos += b.vel*dt}
@body.each{|b| b.vel += b.old_acc*dt}
end
and in rknbody6.rb
def forward(dt)
calc(dt," @old_acc = acc(ba) ")
calc(dt," @pos += @vel*dt ")
calc(dt," @vel += @old_acc*dt ")
end
Bob: The last one is clearly shorter.
Alice: That I don’t mind so much. I’m just not happy with the fact that it is not clear, for a casual reader, what that variable dt is doing there, as the first two arguments of calc, and the appearance of ba is also a mystery; there is no indication here that ba stands for `body array’ and will get its value from @body, somewhere else. In contrast, the earlier version has nothing hidden: the @body.each{|b| . . . } construct is vanilla flavor for someone familiar with Ruby.
Bob: Perhaps we can improve the calc method of Nbody further. How about redefining it in such a way that we can leave out the first argument altogether?
Alice: Ah, that’s a good idea. It is also a logical next step, after introducing the shortcut notion of sending a string in the first place. Once we do something that is somewhat dirty and not so self-explanatory, we might as well go all the way.
Bob: I suppose that we would have to introduce an extra instance variable @dt for Nbody. Otherwise it will not be possible to remove the first argument dt from the current calc. In fact, that would make the definition even shorter. So it would look very clean, like this:
def forward
calc(" @old_acc = acc(ba) ")
calc(" @pos += @vel*dt ")
calc(" @vel += @old_acc*dt ")
end
Note that I have created yet another version of our code, in file rknbody7.rb.
Alice: That is short and sweet, indeed. And we have to modify calc on the Nbody level from what we had before:
def calc(y,s)
@body.each{|b| b.calc(@body,y,s)}
end
to a version with only one parameter, namely the command string:
def calc(s)
@body.each{|b| b.calc(@body, @dt, s)}
end
Since the other two parameters to the calc method of Body are now both instance variables, their value is known here.
Bob: What else do we have to change? We have to set the value of the new variable @dt in the evolve method, and we have to leave out the dt argument when invoking the integration methods. Two small changes, which leaves us with evolve looking like this:
def evolve(integration_method, dt, dt_dia, dt_out, dt_end)
nsteps = 0
e_init
write_diagnostics(nsteps)
@dt = dt
t_dia = dt_dia - 0.5*@dt
t_out = dt_out - 0.5*@dt
t_end = dt_end - 0.5*@dt
while @time < t_end
send(integration_method)
@time += @dt
nsteps += 1
if @time >= t_dia
write_diagnostics(nsteps)
t_dia += dt_dia
end
if @time >= t_out
simple_print
t_out += dt_out
end
end
end
Alice: And the only other changes, besides the change in the forward Euler algorithm we already saw, are the simplified readings of the three other integrators. They now become:
def leapfrog
calc(" @vel += acc(ba)*0.5*dt ")
calc(" @pos += @vel*dt ")
calc(" @vel += acc(ba)*0.5*dt ")
end
def rk2
calc(" @old_pos = @pos ")
calc(" @half_vel = @vel + acc(ba)*0.5*dt ")
calc(" @pos += @vel*0.5*dt ")
calc(" @vel += acc(ba)*dt ")
calc(" @pos = @old_pos + @half_vel*dt ")
end
def rk4
calc(" @old_pos = @pos ")
calc(" @a0 = acc(ba) ")
calc(" @pos = @old_pos + @vel*0.5*dt + @a0*0.125*dt*dt ")
calc(" @a1 = acc(ba) ")
calc(" @pos = @old_pos + @vel*dt + @a1*0.5*dt*dt ")
calc(" @a2 = acc(ba) ")
calc(" @pos = @old_pos + @vel*dt + (@a0+@a1*2)*(1/6.0)*dt*dt ")
calc(" @vel += (@a0+@a1*4+@a2)*(1/6.0)*dt ")
end
Bob: Time to check whether the new code does the same thing as all the older ones:
|gravity> ruby rknbody7a_driver.rb < figure8.in
dt = 0.001
dt_dia = 2.1088
dt_out = 2.1088
dt_end = 2.1088
method = rk4
at time t = 0, after 0 steps :
E_kin = 1.21 , E_pot = -2.5 , E_tot = -1.29
E_tot - E_init = 0
(E_tot - E_init) / E_init = -0
at time t = 2.109, after 2109 steps :
E_kin = 1.21 , E_pot = -2.5 , E_tot = -1.29
E_tot - E_init = -2e-15
(E_tot - E_init) / E_init = 1.55e-15
3
2.1089999999998787e+00
1.0000000000000000e+00
-1.6047303546488470e-04 -1.9320664965417420e-04
-9.3227640249930266e-01 -8.6473492670753516e-01
1.0000000000000000e+00
9.7020367429337440e-01 -2.4296620300772800e-01
4.6595057278750124e-01 4.3244644507801255e-01
1.0000000000000000e+00
-9.7004320125790211e-01 2.4315940965738195e-01
4.6632582971180025e-01 4.3228848162952316e-01
Alice: And so it does. Good! I think we can stick with this version as our work horse, for a while at least.
Bob: There is one improvement I would like to add. Remember when we were debugging our code and we were analyzing a single time step, trying to figure out whether the position and velocity were updated correctly?
Alice: Yes, that was a simple situation, in the case of the forward Euler algorithm, and we could analytically calculate the acceleration.
Bob: Of course, in general we can’t be so lucky, and I would prefer to have an easy way to ask for an output of the acceleration as well. Ah, come to think of it, we have this pretty print output version, that we wrote long ago, but that never used since. Here it is:
def pp # pretty print
print " N = ", @body.size, "\n"
print " time = ", @time, "\n"
@body.each{|b| b.pp}
end
All it does is invoke the pretty print version for each particle:
def pp # pretty print
print to_s
end
which does nothing else put print the Body instance in its standard string form:
def to_s
" mass = " + @mass.to_s + "\n" +
" pos = " + @pos.join(", ") + "\n" +
" vel = " + @vel.join(", ") + "\n"
end
Well, let us just add an extra line for the acceleration. Since pretty print became pp, a pretty print with extras should be called a ppx:
def ppx(body_array) # pretty print, with extra information (acc)
STDERR.print to_s + " acc = " + acc(body_array).join(", ") + "\n"
end
I will put this version in file rknbody8.rb.
Alice: That makes a lot of sense. You choose the standard error output stream STDERR because you don’t want this extra debugging information to be mixed up with the particle output snapshot, which is written on the default standard output stream, and which can be piped into another program.
And of course, you have to give ppx a parameter, namely the array of all the particles in the N-body system, otherwise our particle cannot compute the acceleration that it receives from all other particles. Let’s see. This means that ppx on the Body level must be invoked from the Nbody as:
def ppx # pretty print, with extra information (acc)
print " N = ", @body.size, "\n"
print " time = ", @time, "\n"
@body.each{|b| b.ppx(@body)}
end
Bob: Indeed. Now what else is there left to do? We don’t want to have this extra information all the time, since it would clutter up the output. Let us introduce a special flag, x_flag, a boolean variable that will be set to `true’ if the extra output is desired, and set to `false’ when we don’t need it.
This means that the method that writes the diagnostics will now get x_flag as a second parameter, and a few extra lines at the end, to invoke ppx if the flag is set to be `true’:
def write_diagnostics(nsteps, x_flag)
etot = ekin + epot
STDERR.print <<END
at time t = #{sprintf("%g", time)}, after #{nsteps} steps :
E_kin = #{sprintf("%.3g", ekin)} ,\
E_pot = #{sprintf("%.3g", epot)} ,\
E_tot = #{sprintf("%.3g", etot)}
E_tot - E_init = #{sprintf("%.3g", etot - @e0)}
(E_tot - E_init) / E_init = #{sprintf("%.3g", (etot - @e0)/@e0 )}
END
if x_flag
STDERR.print " for debugging purposes, here is the internal data ",
"representation:\n"
ppx
end
end
Alice: Now that we are adding features, I would like to have the option to echo the initial snapshot, the data file that is read in before any integration step is taken. I know that I can always read that information myself from the input file, but sometimes it is nice to have it right in front of you, together with the other data,
Bob: Yes, especially when are debugging and you want to take a single time step. Okay, let us introduce a second flag init_out. If the value of this flag is `true’, we will require an initial output, in the form of a snapshot on the standard output stream. If the value is `false’, we skip the initial output.
We can implement the effects of both flags, x_flag and init_out, by passing them as additional parameters to evolve. The modifications to the body of this method are then very minor. We only have to change three lines. First, the two calls to write_diagnostics acquire x_flag as extra parameter, as we have already seen. Second, we add an extra line
simple_print if init_out
just before we start the while loop. The new version of the evolve method thus becomes:
def evolve(integration_method, dt, dt_dia, dt_out, dt_end, init_out, x_flag)
nsteps = 0
e_init
write_diagnostics(nsteps, x_flag)
@dt = dt
t_dia = dt_dia - 0.5*dt
t_out = dt_out - 0.5*dt
t_end = dt_end - 0.5*dt
simple_print if init_out
while @time < t_end
send(integration_method)
@time += dt
nsteps += 1
if @time >= t_dia
write_diagnostics(nsteps, x_flag)
t_dia += dt_dia
end
if @time >= t_out
simple_print
t_out += dt_out
end
end
end
Alice: An alternative would have been to make both flags into instance variables for the Nbody class.
Bob: Yes, that would perhaps look more tidy, giving us fewer arguments to pass around. Either way, I don’t care very much: most class definitions involve a trade off between the number of instance variables and the number of arguments being passed around. I’m happy just to keep them as arguments for now.
Alice: Let’s make the necessary changes in the driver file as well:
require "rknbody8.rb" include Math
dt = 0.001 # time step
dt_dia = 2.1088 # diagnostics printing interval
dt_out = 2.1088 # output interval
dt_end = 2.1088 # duration of the integration
init_out = false # initial output requested ?
x_flag = false # extra diagnostics requested ?
#method = "forward" # integration method
#method = "leapfrog" # integration method
#method = "rk2" # integration method
method = "rk4" # integration method
STDERR.print "dt = ", dt, "\n",
"dt_dia = ", dt_dia, "\n",
"dt_out = ", dt_out, "\n",
"dt_end = ", dt_end, "\n",
"init_out = ", init_out, "\n",
"x_flag = ", x_flag, "\n",
"method = ", method, "\n"
nb = Nbody.new nb.simple_read nb.evolve(method, dt, dt_dia, dt_out, dt_end, init_out, x_flag)
Bob: All very straightforward indeed: at the end the two extra parameters for evolve, and above the introduction of the two flags in analogy with the other parameters. Let’s test it out, to see whether we still get the same results for our figure out orbit:
|gravity> ruby rknbody8a_driver.rb < figure8.in
dt = 0.001
dt_dia = 2.1088
dt_out = 2.1088
dt_end = 2.1088
init_out = false
x_flag = false
method = rk4
at time t = 0, after 0 steps :
E_kin = 1.21 , E_pot = -2.5 , E_tot = -1.29
E_tot - E_init = 0
(E_tot - E_init) / E_init = -0
at time t = 2.109, after 2109 steps :
E_kin = 1.21 , E_pot = -2.5 , E_tot = -1.29
E_tot - E_init = -2e-15
(E_tot - E_init) / E_init = 1.55e-15
3
2.1089999999998787e+00
1.0000000000000000e+00
-1.6047303546488470e-04 -1.9320664965417420e-04
-9.3227640249930266e-01 -8.6473492670753516e-01
1.0000000000000000e+00
9.7020367429337440e-01 -2.4296620300772800e-01
4.6595057278750124e-01 4.3244644507801255e-01
1.0000000000000000e+00
-9.7004320125790211e-01 2.4315940965738195e-01
4.6632582971180025e-01 4.3228848162952316e-01
Alice: All is well, clearly. Good! So now we have a version of the code that is both lean and easy to read on the level of the integrators, and has extra debugging options. Progress!.
| Previous | ToC | Up | Next |