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_Force(Body * a, System * s)
  41 {
  42         double distance;
  43         double con;
  44         double gd;
  45
  46         for (unsigned i = 0; i < DIMENSIONS; ++i) //Set Force to zero
  47                 a->F[i] = 0;
  48
  49         for (unsigned index = 0; index < s->N; ++index) // Step through all bodies in the System
  50         {
  51                 Body * b = s->body+index; //Current body
  52                 if (b == a)
  53                         continue; //Otherwise we would have infinite force
  54
  55                 //Calculate distance between a and b
  56                 distance = 0.0;
  57                 for (unsigned i = 0; i < DIMENSIONS; ++i)
  58                         distance += square(b->x[i] - a->x[i]);
  59                 distance = sqrt(distance);
  60                 con = G * a->mass * b->mass / square(distance);
  61                 gd = con / distance;
  62                 for (unsigned i = 0; i < DIMENSIONS; ++i) //Add force from b to a's net force
  63                         a->F[i] += gd * (b->x[i] - a->x[i]);
  64         }
  65 }
  66
  67 /**
  68  * @function Body_Velocity
  69  * @purpose Compute velocity of a body
  70  * @param a - The Body
  71  */
  72 inline void Body_Velocity(Body * a)
  73 {
  74         for (unsigned i = 0; i < DIMENSIONS; ++i)
  75                 a->v[i] += a->F[i] / a->mass * DELTA_T;
  76 }
  77
  78 /**
  79  * @function Body_Position
  80  * @purpose Compute position of a body
  81  * @param a - The Body
  82  */
  83 inline void Body_Position(Body * a)
  84 {
  85         for (unsigned i = 0; i < DIMENSIONS; ++i)
  86                 a->x[i] += a->v[i] * DELTA_T;
  87 }
  88
  89 /**
  90  * @function System_Compute
  91  * @purpose Compute forces and then positions for bodies in System
  92  * NOTE: Only used in the single threaded version of the program
  93  */
  94 inline void System_Compute(System * s)
  95 {
  96         System_Forces(s, s);
  97         System_Positions(s);
  98 }
  99
 100 /**
 101  * @function System_Forces
 102  * @purpose Compute forces on each object in System s1 due to all bodies in System s2
 103  * @param s1, s2 - The two systems
 104  * NOTE: In single threaded version, called with s1 == s2 == &universe
 105  *       In multi threaded version, called with s2 == &universe, but s1 is constructed for each thread
 106  *              In nested multi-threaded version, s2 is constructed for the nested threads as well
 107  */
 108 inline void System_Forces(System * s1, System * s2)
 109 {
 110         //printf("Compute forces for %p - %d bodies...\n", (void*)s1, s1->N);
 111         for (unsigned i = 0; i < s1->N; ++i)
 112         {
 113                 Body_Force(s1->body+i, s2);
 114                 Body_Velocity(s1->body+i);
 115         }
 116         //printf("Compute positions for %p - Done!\n", (void*)s1);
 117 }
 118
 119 /**
 120  * @function System_Positions
 121  * @purpose Compute positions of all objects in System
 122  * @param s - The system
 123  */
 124 inline void System_Positions(System * s)
 125 {
 126         //printf("Compute positions for %p - %d bodies...\n", (void*)s, s->N);
 127         for (unsigned i = 0; i < s->N; ++i)
 128                 Body_Position(s->body+i);
 129         //printf("Compute positions for %p - Done!\n", (void*)s);
 130
 131 }
 132
 133
 134
 135
 136 /**
 137  * @function System_Init
 138  * @purpose Initialise a System from an input file
 139  * @param s - The System
 140  * @param fileName - The input file
 141  */
 142 inline void System_Init(System * s, const char *fileName)
 143 {
 144         char line[LINE_SIZE];
 145         char * token;
 146         FILE * file;
 147
 148         file = fopen(fileName, "rt");
 149         s->N = atoi(fgets(line, LINE_SIZE, file));
 150         s->body = (Body*) calloc((size_t)s->N, sizeof(Body));
 151         s->steps = 0;
 152
 153         //printf("----------------------Initial field-------------------------------\n");
 154
 155         for (unsigned i = 0; i < s->N; ++i)
 156         {
 157                 if (fgets(line, LINE_SIZE, file) != NULL)
 158                 {
 159                         Body * a = s->body+i;
 160                         token = strtok(line, ",; \t");
 161                         a->mass = atof(token);
 162
 163                         for (unsigned j = 0; j < DIMENSIONS; ++j)
 164                         {
 165                                 token = strtok(NULL, ",; \t");
 166                                 a->x[j] = atof(token);
 167                         }
 168                         for (unsigned j = 0; j < DIMENSIONS; ++j)
 169                         {
 170                                 token = strtok(NULL, ",; \t");
 171                                 a->v[j] = atof(token);
 172                         }
 173
 174                         // Ignore any additional tokens
 175                         while (token != NULL)
 176                                 token = strtok(NULL, ",; \t");
 177                         //Body_Print(a);
 178                 }
 179         }
 180
 181         //printf("--------------------------------------------------------------\n");
 182
 183         fclose(file);
 184 }
 185
 186 /**
 187  * @function Universe_Cleanup
 188  * @purpose Cleans up the universe by freeing all memory
 189  *       Also prints the bodies in the universe to a file specified in "options" if required.
 190  */
 191 inline void Universe_Cleanup()
 192 {
 193         //fprintf(stderr, "Called Universe_Cleanup()\n");
 194         if (options.output != NULL) // An output file was specified...
 195         {
 196                 FILE * save = NULL;
 197                 if (strcmp("stdout", options.output) == 0)
 198                         save = stdout;
 199                 else if (strcmp("stderr", options.output) == 0)
 200                         save = stderr;
 201                 else
 202                         save = fopen(options.output, "w");
 203                 if (save == NULL)
 204                         perror("Couldn't save universe to file.");
 205                 else
 206                 {
 207                         // Print the output in the same format as the initial field file
 208                         fprintf(save, "%u\n", universe.N);
 209
 210                         for (unsigned i = 0; i < universe.N; ++i) // Print all bodies to the file
 211                         {
 212                                 Body_Print(universe.body+i, save);
 213
 214
 215                         }
 216                         fclose(save);
 217                 }
 218
 219         }
 220         free(universe.body); // Cleanup memory used for bodies
 221
 222 }
 223
 224 /**
 225  * @function DisplayStatistics
 226  * @purpose Display useful information about the computations done,
 227  *      - Total number of steps computed
 228  *      - Time in seconds elapsed
 229  *      - Clock cycles executed
 230  *      - Equivelant time for a single thread to execute same number of clock cycles
 231  */
 232 inline void DisplayStatistics()
 233 {
 234         clock_t end = clock();
 235         struct timeval end_time;
 236         if (gettimeofday(&end_time, NULL) != 0)
 237         {
 238                 perror("Couldn't get time of day.\n");
 239                 return;
 240         }
 241         float runtime = (float)(end_time.tv_sec-options.start_time.tv_sec) + (float)(end_time.tv_usec-options.start_time.tv_usec) / 1.0e6;
 242         float clock_time = (float)(end - options.start_clock) / (float)(CLOCKS_PER_SEC);
 243         printf("%u\t%f\t%u\t%f\n", universe.steps, runtime, (unsigned)(end - options.start_clock), clock_time);
 244 }
 245
 246
 247 /**
 248  * @function ExitCondition
 249  * @purpose Helper to check whether the program is supposed to exit
 250  *      Does not check for user triggered quit
 251  *      Checks for either a timeout, or computation of the required number of steps
 252  *
 253  * The reason timeouts are integers is because this function is called a lot in lots of threads
 254  *      So using floats (and working out microseconds every time) is expensive
 255  *
 256  */
 257 inline bool ExitCondition(void)
 258 {
 259         bool result = (runstate != RUN ||
 260                         (options.num_steps >= 0 && universe.steps >= options.num_steps) ||
 261                         (options.timeout >= 0 && ((int)(time(NULL) - options.start_time.tv_sec) >= options.timeout)));
 262         //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));
 263         return result;
 264 }