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