Previous | ToC | Up | Next |

** Bob**: Now that we have a shared time step code `nbody_sh1.rb`,
I'd like to make sure that it really gives the same results as the
previous code `nbody_cst1.rb` that used constant time steps.
Having good energy accuracy is one thing, but it would be much better
to check the actual particle positions.

** Alice**: How about computing the distance between the results of two runs,
in 6N-dimensional phase space? A whole N-body system can be viewed
as a single point in a space that spans all the degrees of freedom of
the positions and velocities of all particles.

More generally, if we model a system in D dimensions, a single particle will have D components for its position as well as its velocity, or 2D in total. N particles will have 2DN such degrees of freedom. For planar N-body systems, we need a 4N-dimensional phase space, and for N-body systems in the real three-dimensional world we are dealing with a 6N-dimensional phase space.

** Bob**: In practice this means, I presume, that you use Pythagoras to
compute the distance in such a large space. In other words you take
the square root of the sum of the squares of all the differences, right?

** Alice**: Indeed.

** Bob**: So another way of looking at your suggestion is that you compute
the average distance between corresponding particles in the usual
6-dimensional phase space. Apart from a normalization factor N, the
result would be the same.

** Alice**: Yes, that's right. I like the more abstract picture, but your
description boils down to the same thing.

** Bob**: I find it a lot easier to think about 2D dimensions than about 2ND
dimensions. In either case, the practical procedure is that we have to
subtract two N-body systems, point by point, and then we have to compute
the norm or `size' of the resulting system, a type of absolute value.

** Alice**: Ah, so you * do* like abstract thinking, after all! How about
defining a minus operator, that allows us to subtract two N-body systems,
`nb1` and `nb2`, just by writing `nb1-nb2`? That
would give us a nicely compact notation. And the 2DN-dimensional distance
between the two systems could then be given by `(nb1-nb2).abs`
if we define the correct absolute value operator ` abs` for the `NBody`
class.

** Bob**: Good idea! However, there is no guarantee that two N-body systems
will show their particles in the output in the exact same order. In other
words, what we need is a way to identify which particle in the one system
corresponds to which particle in the other system.

** Alice**: I guess we have to number them, by giving each body a unique
identifier.

** Bob**: Given that every object in Ruby already has an identifier called
`object_id`, a natural name for such an ID would be `body_id`.
Now that is easy to implement. How about writing a nifty tool like this,
in file `nbody_set_id.rb`:

#!/usr/local/bin/ruby -w require "nbody.rb" options_text= <<-END Description: Takes an N-body system, and gives each body a unique ID Long description: This program accepts an N-body system, and assigns a number to each body consecutively, as an instance variable @body_id which takes on integer values. (c) 2005, Piet Hut, Jun Makino; see ACS at www.artcompsi.org example: kali #{$0} -n 0 Short name: -n Long name: --starting_number Value type: int Default value: 1 Description: value of @body_id for 1st body Variable name: n Long description: This option allows the user to start numbering the particles at a value different from the default value (0, say, or 10, or whatever). END c = parse_command_line(options_text) class Body attr_accessor :body_id end nb = ACS_IO.acs_read(NBody) i = c.n nb.body.each do |b| b.body_id = i i += 1 end nb.acs_write($stdout, false, c.precision)

You see, I even allowed for taste: the default is to start numbering at one, but you can start at zero as well, or at any other place you like. Let's test it with our Plummer Model:

|gravity> kali nbody_set_id.rb -h Takes an N-body system, and gives each body a unique ID -n --starting_number : value of @body_id for 1st body [default: 1] --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) |gravity> kali mkplummer.rb -n 2 | kali nbody_set_id.rb --precision 5 ==> Plummer's Model Builder <== Number of particles : N = 2 pseudorandom number seed given : 0 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 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 5 Incremental indentation : add_indent = 2 actual seed used : 252462866881473732163958989257426001914 ACS NBody Array body Body body[0] Fixnum body_id 1 Float mass 5.00000e-01 Vector pos 1.94755e-01 9.73956e-02 -1.22820e-01 Vector vel 1.78002e-01 5.49335e-01 -4.08101e-01 Body body[1] Fixnum body_id 2 Float mass 5.00000e-01 Vector pos -1.94755e-01 -9.73956e-02 1.22820e-01 Vector vel -1.78002e-01 -5.49335e-01 4.08101e-01 Float time 0.00000e+00 SCA |gravity> kali mkplummer.rb -n 2 | kali nbody_set_id.rb --precision 5 -n 42 ==> Plummer's Model Builder <== Number of particles : N = 2 pseudorandom number seed given : 0 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 actual seed used : 210733041374937108452252448370138453664 ==> Takes an N-body system, and gives each body a unique ID <== value of @body_id for 1st body : n = 42 Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 5 Incremental indentation : add_indent = 2 ACS NBody Array body Body body[0] Fixnum body_id 42 Float mass 5.00000e-01 Vector pos 1.16550e-01 -4.85085e-02 2.15785e-01 Vector vel -5.72883e-02 3.45794e-01 6.14121e-01 Body body[1] Fixnum body_id 43 Float mass 5.00000e-01 Vector pos -1.16550e-01 4.85085e-02 -2.15785e-01 Vector vel 5.72883e-02 -3.45794e-01 -6.14121e-01 Float time 0.00000e+00 SCA

** Bob**: Now, with identifiers in place, we can extend our `nbody.rb` file
by adding a subtraction operator:

def -(other) if other.class != self.class raise "other.class.name = #{other.class.name} != #{self.class.name}" end if (n = other.body.size) != body.size raise "other.body.size = #{other.body.size} != #{body.size}" end nb = NBody.new body.each_index do |i| if (id = body[i].body_id) == nil raise "body[#{i}].body_id == nil" end ob = other.body.find{|oi| oi.body_id == id} if ob == nil raise "body_id = #{body[i].body_id} not found in other N-body system" end nb.body[i] = Body.new nb.body[i].pos = body[i].pos - ob.pos nb.body[i].vel = body[i].vel - ob.vel end nb end

This includes proper testing of particle identities.

** Alice**: And I see that you are returning the newly constructed nbody system
` nb` at the end of the method. Also, you're getting quite defensive
in your programming!

** Bob**: This is a situation where defensive programming is called for.
In general, I don't like to litter my code with many error checking
lines all over the place, but in this particular case, it is quite
likely that the difference operator will be invoked for incompatible
N-body systems. So I'm checking whether the other object is really
of ` NBody` class, whether it has the same number of particles as the
calling ` NBody` object, and most importantly, I am then checking
whether each particle in the calling system has a counter part in the
other system, with the same `body_id`.

** Alice**: Applause, applause! As you know, I like to check things, as
a matter of habit, in order to prevent future surprises that may be
hard to debug. But I agree, you can overdo with anything, and there
is a lot to say for keeping code small. In this case, though, the
whole code is so small that it doesn't distract to have these extra
safeguards built in.

** Bob**: With our new tools in hand, it now becomes really trivial to
check for the actual size of the difference between two N-body
systems, by defining the absolute value method, also within the file
`nbody.rb`:

def abs a = 0 body.each{ |b| a += b.pos*b.pos + b.vel*b.vel } sqrt a end

** Alice**: Sometimes we may want to look only at the positions, restricting
ourselves to configuration space, rather than to phase space.

** Bob**: That's easy. Here is a position-only version of ` abs`:

def abs_pos a = 0 body.each{ |b| a += b.pos*b.pos } sqrt a end

** Alice**: Let me try to implement the module that gives the difference
between two N-body systems, as the 6N-dimensional distance in phase space,
or the 2DN-dimensional distance for the general case of D spatial dimensions.
It's time to practice my Ruby skills again. Let's call the file
`nbody_diff.rb`:

#!/usr/local/bin/ruby -w require "nbody.rb" options_text= <<-END Description: 6N-dimensional phase space distance between two N-body systems Long description: This program accepts two N-body systems, and computes the distance between them in the 2DN-dimensional phase space of all positions and velocities, if they are given in D spatial dimensions. The values of the masses of the particles are not used, and each particle is given equal weight. In order to compare corresponding particles in the two systems, the body_id instance variable of each body is used as an identifier. (c) 2005, Piet Hut, Jun Makino; see ACS at www.artcompsi.org example: cat run1.out run2.out | kali #{$0} Short name: -r Long name: --positions_only Value type: bool Description: Using only positions, not velocities Variable name: position_only_flag Long description: This option allows the user to compute the distance between two N-body systems using only particle information; velocity information is not used. END c = parse_command_line(options_text) class Body attr_accessor :body_id end nb1 = ACS_IO.acs_read(NBody) nb2 = ACS_IO.acs_read(NBody) nb = nb1 - nb2 d = nb.body[0].pos.size n = nb.body.size print "#{2*d}N-dim. phase space dist. for two #{n}-body systems: " if c.position_only_flag == false printf(" %.#{c.precision}e\n", nb.abs) else printf(" %.#{c.precision}e\n", nb.abs_pos) end

How about this?

** Bob**: Looks fine, as far as I can see. You put in options to allow for more
precision, without burdening the reader with all those extra digits in
the general case. And you also gave an option for positions-only
comparisons.

Time to test it! We still have the initial conditions lying around
for a circular binary, I believe. Ah yes, here it is, in
`circular_binary.in`:

ACS NBody Array body Body body[0] Float mass 4 Vector pos 1 0 Vector vel 0 1 Body body[1] Float mass 4 Vector pos -1 0 Vector vel 0 -1 SCATo start with, let's see whether two identical systems are assigned a distance of zero correctly:

|gravity> kali nbody_set_id.rb < circular_binary.in > tmp0.in ==> Takes an N-body system, and gives each body a unique ID <== value of @body_id for 1st body : 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 |gravity> cat tmp0.in tmp0.in | 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 4N-dim. phase space dist. for two 2-body systems: 0.0000000000000000e+00

ACS NBody Array body Body body[0] Float mass 4 Vector pos 1 0.03 Vector vel 0 1 Body body[1] Float mass 4 Vector pos -1 0 Vector vel 0 -1.04 SCAAnd here is our first real result:

|gravity> kali nbody_set_id.rb < pert_circ_binary.in > tmp1.in ==> Takes an N-body system, and gives each body a unique ID <== value of @body_id for 1st body : 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 |gravity> cat tmp0.in tmp1.in | 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 4N-dim. phase space dist. for two 2-body systems: 5.0000000000000031e-02Just what it should be! Now for positions only, and a higher precision:

|gravity> cat tmp0.in tmp1.in | kali nbody_diff.rb -r --precision 5 ==> 6N-dimensional phase space distance between two N-body systems <== Using only positions, not velocities Screen Output Verbosity Level : verbosity = 1 ACS Output Verbosity Level : acs_verbosity = 1 Floating point precision : precision = 5 Incremental indentation : add_indent = 2 4N-dim. phase space dist. for two 2-body systems: 3.00000e-02Perfect. I think we're ready to move on.

Previous | ToC | Up | Next |