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         float distance;
  66         float con;
  67         float gd;
  68         //Calculate distance between a and b
  69         distance = 0.0;
  70         for (unsigned i = 0; i < DIMENSIONS; ++i)
  71                 distance += square(b->x[i] - a->x[i]);
  72         distance = sqrt(distance);
  73         con = G * a->mass * b->mass / square(distance);
  74         gd = con / distance;
  75         for (unsigned i = 0; i < DIMENSIONS; ++i) //Add force from b to a's net force
  76                 a->F[i] += gd * (b->x[i] - a->x[i]);
  77 }
  78
  79 /**
  80  * @function Body_Velocity
  81  * @purpose Compute velocity of a body
  82  * @param a - The Body
  83  */
  84 inline void Body_Velocity(Body * a)
  85 {
  86         for (unsigned i = 0; i < DIMENSIONS; ++i)
  87                 a->v[i] += a->F[i] / a->mass * DELTA_T;
  88 }
  89
  90 /**
  91  * @function Body_Position
  92  * @purpose Compute position of a body
  93  * @param a - The Body
  94  */
  95 inline void Body_Position(Body * a)
  96 {
  97         for (unsigned i = 0; i < DIMENSIONS; ++i)
  98                 a->x[i] += a->v[i] * DELTA_T;
  99 }
 100
 101 /**
 102  * @function System_Compute
 103  * @purpose Compute forces and then positions for bodies in System
 104  * NOTE: Only used in the single threaded version of the program
 105  */
 106 inline void System_Compute(System * s)
 107 {
 108         System_Forces(s, s);
 109         System_Positions(s);
 110 }
 111
 112 /**
 113  * @function System_Forces
 114  * @purpose Compute forces on each object in System s1 due to all bodies in System s2
 115  * @param s1, s2 - The two systems
 116  * NOTE: In single threaded version, called with s1 == s2 == &universe
 117  *       In multi threaded version, called with s2 == &universe, but s1 is constructed for each thread
 118  *              In nested multi-threaded version, s2 is constructed for the nested threads as well
 119  */
 120 inline void System_Forces(System * s1, System * s2)
 121 {
 122         //printf("Compute forces for %p - %d bodies...\n", (void*)s1, s1->N);
 123         for (unsigned i = 0; i < s1->N; ++i)
 124         {
 125                 Body_Forces(s1->body+i, s2);
 126                 Body_Velocity(s1->body+i);
 127         }
 128         //printf("Compute positions for %p - Done!\n", (void*)s1);
 129 }
 130
 131 /**
 132  * @function System_Positions
 133  * @purpose Compute positions of all objects in System
 134  * @param s - The system
 135  */
 136 inline void System_Positions(System * s)
 137 {
 138         //printf("Compute positions for %p - %d bodies...\n", (void*)s, s->N);
 139         for (unsigned i = 0; i < s->N; ++i)
 140                 Body_Position(s->body+i);
 141         //printf("Compute positions for %p - Done!\n", (void*)s);
 142
 143 }
 144
 145
 146
 147
 148 /**
 149  * @function System_Init
 150  * @purpose Initialise a System from an input file
 151  * @param s - The System
 152  * @param fileName - The input file
 153  */
 154 inline void System_Init(System * s, const char *fileName)
 155 {
 156         char line[LINE_SIZE];
 157         char * token;
 158         FILE * file;
 159
 160         file = fopen(fileName, "rt");
 161
 162         if (file == NULL)
 163         {
 164                 fprintf(stderr, "%s\n", fileName);
 165                 perror("Couldn't open file");
 166                 exit(EXIT_FAILURE);
 167         }
 168
 169         s->N = atoi(fgets(line, LINE_SIZE, file));
 170         s->body = (Body*) calloc((size_t)s->N, sizeof(Body));
 171         s->steps = 0;
 172
 173         //printf("----------------------Initial field-------------------------------\n");
 174
 175         for (unsigned i = 0; i < s->N; ++i)
 176         {
 177                 if (fgets(line, LINE_SIZE, file) != NULL)
 178                 {
 179                         Body * a = s->body+i;
 180                         token = strtok(line, ",; \t");
 181                         a->mass = atof(token);
 182
 183                         for (unsigned j = 0; j < DIMENSIONS; ++j)
 184                         {
 185                                 token = strtok(NULL, ",; \t");
 186                                 a->x[j] = atof(token);
 187                         }
 188                         for (unsigned j = 0; j < DIMENSIONS; ++j)
 189                         {
 190                                 token = strtok(NULL, ",; \t");
 191                                 a->v[j] = atof(token);
 192                         }
 193
 194                         // Ignore any additional tokens
 195                         while (token != NULL)
 196                                 token = strtok(NULL, ",; \t");
 197                         //Body_Print(a);
 198                 }
 199         }
 200
 201         //printf("--------------------------------------------------------------\n");
 202
 203         fclose(file);
 204 }
 205
 206 /**
 207  * @function Universe_Cleanup
 208  * @purpose Cleans up the universe by freeing all memory
 209  *       Also prints the bodies in the universe to a file specified in "options" if required.
 210  */
 211 inline void Universe_Cleanup()
 212 {
 213         //fprintf(stderr, "Called Universe_Cleanup()\n");
 214         if (options.output != NULL) // An output file was specified...
 215         {
 216                 FILE * save = NULL;
 217                 if (strcmp("stdout", options.output) == 0)
 218                         save = stdout;
 219                 else if (strcmp("stderr", options.output) == 0)
 220                         save = stderr;
 221                 else
 222                         save = fopen(options.output, "w");
 223                 if (save == NULL)
 224                         perror("Couldn't save universe to file.");
 225                 else
 226                 {
 227                         // Print the output in the same format as the initial field file
 228                         fprintf(save, "%u\n", universe.N);
 229
 230                         for (unsigned i = 0; i < universe.N; ++i) // Print all bodies to the file
 231                         {
 232                                 Body_Print(universe.body+i, save);
 233
 234
 235                         }
 236                         fclose(save);
 237                 }
 238
 239         }
 240         free(universe.body); // Cleanup memory used for bodies
 241
 242 }
 243
 244 /**
 245  * @function DisplayStatistics
 246  * @purpose Display useful information about the computations done,
 247  *      - Total number of steps computed
 248  *      - Time in seconds elapsed
 249  *      - Clock cycles executed
 250  *      - Equivelant time for a single thread to execute same number of clock cycles
 251  */
 252 inline void DisplayStatistics()
 253 {
 254         clock_t end = clock();
 255         struct timeval end_time;
 256         if (gettimeofday(&end_time, NULL) != 0)
 257         {
 258                 perror("Couldn't get time of day.\n");
 259                 return;
 260         }
 261         float runtime = (float)(end_time.tv_sec-options.start_time.tv_sec) + (float)(end_time.tv_usec-options.start_time.tv_usec) / 1.0e6;
 262         float clock_time = (float)(end - options.start_clock) / (float)(CLOCKS_PER_SEC);
 263         printf("%u\t%f\t%u\t%f\n", universe.steps, runtime, (unsigned)(end - options.start_clock), clock_time);
 264 }
 265
 266
 267 /**
 268  * @function ExitCondition
 269  * @purpose Helper to check whether the program is supposed to exit
 270  *      Does not check for user triggered quit
 271  *      Checks for either a timeout, or computation of the required number of steps
 272  *
 273  * The reason timeouts are integers is because this function is called a lot in lots of threads
 274  *      So using floats (and working out microseconds every time) is expensive
 275  *
 276  */
 277 inline bool ExitCondition(void)
 278 {
 279         bool result = (runstate != RUN ||
 280                         (options.num_steps >= 0 && universe.steps >= options.num_steps) ||
 281                         (options.timeout >= 0 && ((int)(time(NULL) - options.start_time.tv_sec) >= options.timeout)));
 282         //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));
 283         return result;
 284 }
 285
 286 /**
 287  * @function Split_System
 288  * @purpose Helper to divide one system into an array of systems
 289  *      Each sub system will have N = (s->N / n) bodies in it
 290  * @param s - The original system (typically &universe)
 291  * @param n - The number of sub systems in the array
 292  *
 293  * WARNING: It is the caller's responsibility to free() the returned array
 294  */
 295 System * Split_System(System * s, unsigned n)
 296 {
 297         System * result = (System*)(calloc(n, sizeof(System)));
 298         if (result == NULL)
 299         {
 300                 perror("Couldn't create array of sub systems");
 301                 exit(EXIT_FAILURE);
 302         }
 303
 304         unsigned n_per_system = (s->N) / n;
 305         unsigned remainder = (s->N) % n;
 306
 307         for (unsigned i = 0; i < n; ++i)
 308         {
 309                 result[i].N = n_per_system;
 310                 if (i == n-1)
 311                         result[i].N += remainder;
 312                 result[i].body = (s->body) + (n_per_system * i);
 313                 result[i].steps = 0;
 314         }
 315         return result;
 316 }
 317
 318 /**
 319  * @function System_Random
 320  * @purpose Randomly generate initial body field
 321  * @param s - The system
 322  * @param n - Number of bodies
 323  */
 324 void System_Random(System * s, int n)
 325 {
 326         srand(time(NULL));
 327         s->N = (unsigned)n;
 328         s->body = (Body*) calloc((size_t)s->N, sizeof(Body));
 329         s->steps = 0;
 330
 331         s->body[0].mass = 1e11;
 332
 333         for (unsigned a = 1; a < s->N; ++a)
 334         {
 335                 s->body[a].mass = 1e7;
 336                 for (unsigned i = 0; i < DIMENSIONS; ++i)
 337                 {
 338                         s->body[a].x[i] = -10000 + rand() % 20000;
 339                         s->body[a].v[i] = -2 + rand() % 4;
 340                 }
 341                 s->body[a].v[2] = 0;
 342                 s->body[a].x[2] = 0; //set z to zero; force in plane
 343
 344         }
 345 }