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