course/semester2/pprog/assignment1/single-thread/nbody.c~

   1 /**
   2  * @file nbody.c
   3  * @author Sam Moore (20503628) 2012
   4  * @purpose N-Body simulator - Definition of simulation functions; single threaded version
   5  */
   6
   7 #include <stdlib.h>
   8 #include <stdio.h>
   9 #include <math.h>
  10 #include <time.h>
  11 #include <string.h>
  12 #include <stdbool.h>
  13 #include <stdlib.h>
  14 #include <unistd.h>
  15
  16 #include "nbody.h"
  17
  18 RUNSTATE runstate = RUN;
  19
  20 /**
  21  * @function Body_Print
  22  * @purpose Print the body to a file (or stdout)
  23  * @param a - Body to print'
  24  * @param out - FILE* to print to
  25  */
  26 inline void Body_Print(Body * a, FILE * out)
  27 {
  28         //fprintf(out,"Body %p M=%f X=%f Y=%f Z=%f Fx=%f Fy=%f Fz=%f Vx=%f Vy=%f Vz=%f\n",
  29            //(void*)a, a->mass, a->x[0], a->x[1], a->x[2], a->F[0], a->F[1], a->F[2], a->v[0], a->v[1], a->v[2]);
  30         fprintf(out, "%f %f %f %f %f %f %f %f %f %f\n",
  31                 a->mass, a->x[0], a->x[1], a->x[2], a->v[0], a->v[1], a->v[2], a->F[0], a->F[1], a->F[2]);
  32 }
  33
  34 /**
  35  * @function Body_Force
  36  * @purpose Calculates the force on a single Body due to all bodies in a System
  37  * @param a - The Body
  38  * @param b - The System
  39  */
  40 inline void Body_Forces(Body * a, System * s)
  41 {
  42
  43
  44         for (unsigned i = 0; i < DIMENSIONS; ++i) //Set Force to zero
  45                 a->F[i] = 0;
  46
  47         for (unsigned index = 0; index < s->N; ++index) // Step through all bodies in the System
  48         {
  49                 Body * b = s->body+index; //Current body
  50                 if (b == a)
  51                         continue; //Otherwise we would have infinite force
  52
  53                 Body_Force(a, b);
  54         }
  55 }
  56
  57 /**
  58  * @function Body_Force
  59  * @purpose Compute force on body a due to body b
  60  * @param a - Body a
  61  * @param b - Body b
  62  */
  63 inline void Body_Force(Body * a, Body * b)
  64 {
  65         //Calculate distance between a and b
  66         distance = 0.0;
  67         for (unsigned i = 0; i < DIMENSIONS; ++i)
  68                 distance += square(b->x[i] - a->x[i]);
  69         distance = sqrt(distance);
  70         con = G * a->mass * b->mass / square(distance);
  71         gd = con / distance;
  72         for (unsigned i = 0; i < DIMENSIONS; ++i) //Add force from b to a's net force
  73                 a->F[i] += gd * (b->x[i] - a->x[i]);
  74 }
  75
  76 /**
  77  * @function Body_Velocity
  78  * @purpose Compute velocity of a body
  79  * @param a - The Body
  80  */
  81 inline void Body_Velocity(Body * a)
  82 {
  83         for (unsigned i = 0; i < DIMENSIONS; ++i)
  84                 a->v[i] += a->F[i] / a->mass * DELTA_T;
  85 }
  86
  87 /**
  88  * @function Body_Position
  89  * @purpose Compute position of a body
  90  * @param a - The Body
  91  */
  92 inline void Body_Position(Body * a)
  93 {
  94         for (unsigned i = 0; i < DIMENSIONS; ++i)
  95                 a->x[i] += a->v[i] * DELTA_T;
  96 }
  97
  98 /**
  99  * @function System_Compute
 100  * @purpose Compute forces and then positions for bodies in System
 101  * NOTE: Only used in the single threaded version of the program
 102  */
 103 inline void System_Compute(System * s)
 104 {
 105         System_Forces(s, s);
 106         System_Positions(s);
 107 }
 108
 109 /**
 110  * @function System_Forces
 111  * @purpose Compute forces on each object in System s1 due to all bodies in System s2
 112  * @param s1, s2 - The two systems
 113  * NOTE: In single threaded version, called with s1 == s2 == &universe
 114  *       In multi threaded version, called with s2 == &universe, but s1 is constructed for each thread
 115  *              In nested multi-threaded version, s2 is constructed for the nested threads as well
 116  */
 117 inline void System_Forces(System * s1, System * s2)
 118 {
 119         //printf("Compute forces for %p - %d bodies...\n", (void*)s1, s1->N);
 120         for (unsigned i = 0; i < s1->N; ++i)
 121         {
 122                 Body_Forces(s1->body+i, s2);
 123                 Body_Velocity(s1->body+i);
 124         }
 125         //printf("Compute positions for %p - Done!\n", (void*)s1);
 126 }
 127
 128 /**
 129  * @function System_Positions
 130  * @purpose Compute positions of all objects in System
 131  * @param s - The system
 132  */
 133 inline void System_Positions(System * s)
 134 {
 135         //printf("Compute positions for %p - %d bodies...\n", (void*)s, s->N);
 136         for (unsigned i = 0; i < s->N; ++i)
 137                 Body_Position(s->body+i);
 138         //printf("Compute positions for %p - Done!\n", (void*)s);
 139
 140 }
 141
 142
 143
 144
 145 /**
 146  * @function System_Init
 147  * @purpose Initialise a System from an input file
 148  * @param s - The System
 149  * @param fileName - The input file
 150  */
 151 inline void System_Init(System * s, const char *fileName)
 152 {
 153         char line[LINE_SIZE];
 154         char * token;
 155         FILE * file;
 156
 157         file = fopen(fileName, "rt");
 158         s->N = atoi(fgets(line, LINE_SIZE, file));
 159         s->body = (Body*) calloc((size_t)s->N, sizeof(Body));
 160         s->steps = 0;
 161
 162         //printf("----------------------Initial field-------------------------------\n");
 163
 164         for (unsigned i = 0; i < s->N; ++i)
 165         {
 166                 if (fgets(line, LINE_SIZE, file) != NULL)
 167                 {
 168                         Body * a = s->body+i;
 169                         token = strtok(line, ",; \t");
 170                         a->mass = atof(token);
 171
 172                         for (unsigned j = 0; j < DIMENSIONS; ++j)
 173                         {
 174                                 token = strtok(NULL, ",; \t");
 175                                 a->x[j] = atof(token);
 176                         }
 177                         for (unsigned j = 0; j < DIMENSIONS; ++j)
 178                         {
 179                                 token = strtok(NULL, ",; \t");
 180                                 a->v[j] = atof(token);
 181                         }
 182
 183                         // Ignore any additional tokens
 184                         while (token != NULL)
 185                                 token = strtok(NULL, ",; \t");
 186                         //Body_Print(a);
 187                 }
 188         }
 189
 190         //printf("--------------------------------------------------------------\n");
 191
 192         fclose(file);
 193 }
 194
 195 /**
 196  * @function Universe_Cleanup
 197  * @purpose Cleans up the universe by freeing all memory
 198  *       Also prints the bodies in the universe to a file specified in "options" if required.
 199  */
 200 inline void Universe_Cleanup()
 201 {
 202         //fprintf(stderr, "Called Universe_Cleanup()\n");
 203         if (options.output != NULL) // An output file was specified...
 204         {
 205                 FILE * save = NULL;
 206                 if (strcmp("stdout", options.output) == 0)
 207                         save = stdout;
 208                 else if (strcmp("stderr", options.output) == 0)
 209                         save = stderr;
 210                 else
 211                         save = fopen(options.output, "w");
 212                 if (save == NULL)
 213                         perror("Couldn't save universe to file.");
 214                 else
 215                 {
 216                         // Print the output in the same format as the initial field file
 217                         fprintf(save, "%u\n", universe.N);
 218
 219                         for (unsigned i = 0; i < universe.N; ++i) // Print all bodies to the file
 220                         {
 221                                 Body_Print(universe.body+i, save);
 222
 223
 224                         }
 225                         fclose(save);
 226                 }
 227
 228         }
 229         free(universe.body); // Cleanup memory used for bodies
 230
 231 }
 232
 233 /**
 234  * @function DisplayStatistics
 235  * @purpose Display useful information about the computations done,
 236  *      - Total number of steps computed
 237  *      - Time in seconds elapsed
 238  *      - Clock cycles executed
 239  *      - Equivelant time for a single thread to execute same number of clock cycles
 240  */
 241 inline void DisplayStatistics()
 242 {
 243         clock_t end = clock();
 244         struct timeval end_time;
 245         if (gettimeofday(&end_time, NULL) != 0)
 246         {
 247                 perror("Couldn't get time of day.\n");
 248                 return;
 249         }
 250         float runtime = (float)(end_time.tv_sec-options.start_time.tv_sec) + (float)(end_time.tv_usec-options.start_time.tv_usec) / 1.0e6;
 251         float clock_time = (float)(end - options.start_clock) / (float)(CLOCKS_PER_SEC);
 252         printf("%u\t%f\t%u\t%f\n", universe.steps, runtime, (unsigned)(end - options.start_clock), clock_time);
 253 }
 254
 255
 256 /**
 257  * @function ExitCondition
 258  * @purpose Helper to check whether the program is supposed to exit
 259  *      Does not check for user triggered quit
 260  *      Checks for either a timeout, or computation of the required number of steps
 261  *
 262  * The reason timeouts are integers is because this function is called a lot in lots of threads
 263  *      So using floats (and working out microseconds every time) is expensive
 264  *
 265  */
 266 inline bool ExitCondition(void)
 267 {
 268         bool result = (runstate != RUN ||
 269                         (options.num_steps >= 0 && universe.steps >= options.num_steps) ||
 270                         (options.timeout >= 0 && ((int)(time(NULL) - options.start_time.tv_sec) >= options.timeout)));
 271         //fprintf(stderr,"runstate %d\n timeout %d\n steps %d\n", (int)(runstate != RUN), (int)(options.timeout > 0.00 && ((unsigned)(time(NULL) - options.start_time.tv_sec) >= options.timeout)), (int)(options.num_steps > 0 && universe.steps > options.num_steps));
 272         return result;
 273 }
 274
 275 /**
 276  * @function Split_System
 277  * @purpose Helper to divide one system into an array of systems
 278  *      Each sub system will have N = (s->N / n) bodies in it
 279  * @param s - The original system (typically &universe)
 280  * @param n - The number of sub systems in the array
 281  *
 282  * WARNING: It is the caller's responsibility to free() the returned array
 283  */
 284 System * Split_System(System * s, unsigned n)
 285 {
 286         System * result = (System*)(calloc(n, sizeof(System)));
 287         if (result == NULL)
 288         {
 289                 perror("Couldn't create array of sub systems");
 290                 exit(EXIT_FAILURE);
 291         }
 292
 293         unsigned n_per_system = (s->N) / n;
 294         unsigned remainder = (s->N) % n;
 295
 296         for (unsigned i = 0; i < n; ++i)
 297         {
 298                 result[i].N = n_per_system;
 299                 if (i == n-1)
 300                         result[i].N += remainder;
 301                 result[i].body = (s->body) + (n_per_system * i);
 302                 result[i].steps = 0;
 303         }
 304         return result;
 305 }