NBabel

 1 #!/usr/bin/perl
 2 
 3 =head1 NAME
 4 
 5 nbabel - Routine to evolve a simple N-body system using leap-frog integration.
 6 
 7 =head1 SYNOPSIS
 8 
 9 This consists of a number of routines needed to evolve an N-body system with a
10 direct force calculation using a leap-frog integrator. The routine that does the
11 actual
12 evolution is C<evolve_cluster> below which is called with the name of the file
13 containing the initial conditions. Thus a synopsis would be:
14 
15    evolve_cluster('input16')
16 
17 =head1 DESCRIPTION
18 
19 The code reads in initial conditions from a file exected to have an ID column,
20 mass, 3D positions and 3D velocities. These are so evolved forward with a simple
21 leap-frog scheme. See the description of C<evolve_cluster> below.
22 
23 I have not used a fully vectorised approach - rather I have kept the x, y and z
24 components of the positions and velocities in a Perl array. This makes it easier
25 to write the code for the relative separations without going into dimension
26 fiddling and it makes the code more readable. There is a noticeable hit in
27 performance at low N because the (slow) Perl loops over the dimensions is
28 significant. At large N this extra overhead becomes pretty much irrelevant.
29 
30 =cut
31 
32 use strict;
33 use warnings;
34 use PDL;
35 use PDL::NiceSlice;
36 
37 use Time::HiRes qw(usleep ualarm gettimeofday tv_interval);
38 
39 
40 #
41 # Simple N-body calculator in PDL
42 #
43 
44 
45 # Get the initial conditions from the command line
46 my $DIR = '../InputFiles/';
47 my @files = ('input16', 'input32', 'input64', 'input128', 'input256', 'input512')
48             ;
49 @files = ('input1k');
50 
51 
52 foreach my $f (@files) {
53   print "Doing = $f - ";
54   evolve_cluster($DIR.$f);
55 }
56 
57 
58 =head2 evolve_cluster
59 
60 Core routine to evolve the cluster. The procedure is first to read in the
61 initial conditions with C<read_initial_conditions>. The subsequent steps are
62 then
63 
64 =over 4
65 
66 =item Calculate the initial acceleration and energy of the system.
67 
68 This is stored in Store this in C<$a0> and C<$E0> respectively.
69 
70 =item Update the positions
71 
72 =item Calculate the acceleration with this new position but old velocities.
73 
74 =item Update velocities
75 
76 =item Assign current acceleration to C<$a0> and repeat.
77 
78 =item Repeat until t=1 in N-body units
79 
80 =back
81 
82 =cut
83 
84 sub evolve_cluster {
85   my $file = shift;
86 
87   my ($m, $pos, $vel) = read_initial_conditions($file);
88 
89   my $dt = 0.001;
90   my ($a0, $norm) = calculate_acceleration($m, $pos);
91   my $E0 = calculate_energy($m, $pos, $vel, $norm);
92 
93   my $tend = 1.0;
94   my ($nsteps, $tcurrent) = (0, 0.0);
95 
96   my $t1 = [gettimeofday];
97   while ($tcurrent < $tend) {
98     $pos = update_position($pos, $vel, $a0, $dt);
99     ($a, $norm) = calculate_acceleration($m, $pos);
100     $vel = update_velocity($vel, $a0, $a, $dt);
101 
102     $a0 = $a;
103 
104     $tcurrent += $dt;
105     $nsteps++;
106 
107 
108     # This keeps track of the evolution of the system - you could
109     # comment this out for speed when you have very few particles but
110     # it should not give a noticeable overhead for large N.
111     if (($nsteps % 10) == 0) {
112       my $E = calculate_energy($m, $pos, $vel, $norm);
113       printf "t=%.4f Etot=%8.5f  dEtot=%8.5f\n", $tcurrent, $E->at(0), $E->at(0)-
114                 $E0->at(0)#
115    }
116   }
117   my $t2 = [gettimeofday];
118 
119 
120   my $E = calculate_energy($m, $pos, $vel, $norm);
121   printf "Execution time = %7.4fs  - ", tv_interval($t1, $t2);
122   printf "t=%.4f Etot=%8.5f  dEtot=%10.4e\n", $tcurrent, $E->at(0), $E->at(0)-
123             $E0->at(0);
124 
125 }
126 
127 
128 
129 #
130 # Utility functions below
131 
132 #
133 sub calculate_energy {
134   my ($m, $pos, $vel, $norm) = @_;
135 
136   my $E_kin=0.0;
137   for (my $i=0; $i<3; $i++) {
138     $E_kin += $vel->[$i]**2
139   }
140   $E_kin = sum(0.5*$m*$E_kin);
141 
142   my $m_ij = outer($m, $m);
143   my $E_pot = $m_ij/$norm;
144 
145   $E_pot->diagonal(0, 1) .= 0.0; # No self-energy
146   $E_pot = -0.5*sum($E_pot);
147 
148   return pdl($E_kin+$E_pot, $E_kin, $E_pot);
149 
150 }
151 
152 sub calculate_acceleration {
153   my ($m, $pos) = @_;
154 
155   my $n_obj = $m->nelem();
156 
157   # Loop over dimensions
158   my @diff=();
159   my $norm = 0.;
160   for (my $i=0; $i<3; $i++) {
161     $diff[$i] = $pos->[$i]-transpose($pos->[$i]);
162     $diff[$i]->diagonal(0, 1) .= 0.0; # Ignore diagonals
163     $norm += $diff[$i]*$diff[$i];
164   }
165 
166   $norm = sqrt($norm);
167   $norm->diagonal(0, 1) .= 1e30;
168   my $norm3 = $norm**3;
169 
170   my $a = [];
171   for (my $i=0; $i<3; $i++) {
172     my $tmp = -$m*$diff[$i]/$norm3;
173     $tmp->diagonal(0, 1) .= 0;
174     $a->[$i] = sumover($tmp->xchg(0, 1));
175   }
176   return ($a, $norm);
177 
178 }
179 
180 sub update_position {
181   my ($pos, $vel, $a, $dt) = @_;
182   # Loop over dimensions
183   for (my $i=0; $i<3; $i++) {
184     $pos->[$i] += $vel->[$i]*$dt + 0.5*$a->[$i]*$dt*$dt;
185   }
186 
187   return $pos;
188 }
189 
190 sub update_velocity {
191   my ($vel, $a, $a0, $dt) = @_;
192   # Loop over dimensions
193   for (my $i=0; $i<3; $i++) {
194     $vel->[$i] += 0.5*($a->[$i] + $a0->[$i])*$dt;
195   }
196   return $vel;
197 }
198 
199 sub read_initial_conditions {
200   my $file = shift;
201   my ($id, $m, $x, $y, $z, $vx, $vy, $vz) = rcols $file;
202 
203   # Drop empty lines....
204   my $use = which($m > 0);
205   $m = $m($use;|);
206   $x = $x($use;|);
207   $y = $y($use;|);
208   $z = $z($use;|);
209   $vx = $vx($use;|);
210   $vy = $vy($use;|);
211   $vz = $vz($use;|);
212 
213   # Note we use Perl arrays  - this probably could have been handled better
214   return ($m, [$x, $y, $z], [$vx, $vy, $vz]);
215 }
216