course/semester2/pprog/assignment1/mthread/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 //Stuff to do with worker threads
  19 // The worker threads (if used) are spawned by the computation thread, not the main thread
  20 // The main thread (which handles initialisation, graphics, user input, cleanup and exit) does not need to know about these.
  21 #if NUM_THREADS > 1
  22
  23         pthread_t worker_thread[NUM_THREADS]; //Worker threads responsible for Force and Position updates
  24         System sub_system[NUM_THREADS]; //Array of System's used to divide up the main "universe" System for worker threads
  25         pthread_mutex_t mutex_workers;
  26         pthread_cond_t workers_done_cv;
  27         unsigned worker_count;
  28 #endif
  29
  30 /**
  31  * Prints the body on screen
  32  */
  33 void Body_Print(Body * a)
  34 {
  35         printf("Body %p M=%f X=%f Y=%f Z=%f Fx=%f Fy=%f Fz=%f Vx=%f Vy=%f Vz=%f\n",
  36            (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]);
  37 }
  38
  39 /**
  40  * Computing forces
  41  */
  42 void Body_Force(Body * a, System * s)
  43 {
  44         double distance;
  45         double con;
  46         double gd;
  47
  48         for (unsigned i = 0; i < DIMENSIONS; ++i)
  49                 a->F[i] = 0;
  50
  51         for (unsigned index = 0; index < s->N; ++index)
  52         {
  53                 Body * b = s->body+index;
  54                 if (b == a)
  55                         continue;
  56
  57                 distance = 0.0;
  58                 for (unsigned i = 0; i < DIMENSIONS; ++i)
  59                         distance += square(b->x[i] - a->x[i]);
  60                 distance = sqrt(distance);
  61                 con = G * a->mass * b->mass / square(distance);
  62                 gd = con / distance;
  63                 for (unsigned i = 0; i < DIMENSIONS; ++i)
  64                         a->F[i] += gd * (b->x[i] - a->x[i]);
  65         }
  66 }
  67
  68 /**
  69  * Compute velocities
  70  */
  71 void Body_Velocity(Body * a)
  72 {
  73         for (unsigned i = 0; i < DIMENSIONS; ++i)
  74                 a->v[i] += a->F[i] / a->mass * DELTA_T;
  75 }
  76
  77 /**
  78  * Compute positions
  79  */
  80 void Body_Position(Body * a)
  81 {
  82         for (unsigned i = 0; i < DIMENSIONS; ++i)
  83                 a->x[i] += a->v[i] * DELTA_T;
  84 }
  85
  86 /**
  87  * Compute forces on each object in System s1 due to all bodies in System s2
  88  */
  89 void System_Forces(System * s1, System * s2)
  90 {
  91         //printf("Compute forces for %p - %d bodies...\n", (void*)s1, s1->N);
  92         for (unsigned i = 0; i < s1->N; ++i)
  93         {
  94                 Body_Force(s1->body+i, s2);
  95                 Body_Velocity(s1->body+i);
  96         }
  97         //printf("Compute positions for %p - Done!\n", (void*)s1);
  98 }
  99
 100 /**
 101  * Compute positions of all objects in System
 102  */
 103 void System_Positions(System * s)
 104 {
 105         //printf("Compute positions for %p - %d bodies...\n", (void*)s, s->N);
 106         for (unsigned i = 0; i < s->N; ++i)
 107                 Body_Position(s->body+i);
 108         //printf("Compute positions for %p - Done!\n", (void*)s);
 109 }
 110
 111
 112
 113
 114 /*
 115  * This function reads an input file. You can change it if you choose a
 116  * different file format
 117  */
 118 #define LINE_SIZE BUFSIZ
 119 void System_Init(System * s, char *fileName)
 120 {
 121         char line[LINE_SIZE];
 122         char * token;
 123         FILE * file;
 124
 125         file = fopen(fileName, "rt");
 126         s->N = atoi(fgets(line, LINE_SIZE, file));
 127         s->body = (Body*) calloc((size_t)s->N, sizeof(Body));
 128
 129         //printf("----------------------Initial field-------------------------------\n");
 130
 131         for (unsigned i = 0; i < s->N; ++i)
 132         {
 133                 if (fgets(line, LINE_SIZE, file) != NULL)
 134                 {
 135                         Body * a = s->body+i;
 136                         token = strtok(line, ",; ");
 137                         a->mass = atof(token);
 138
 139                         for (unsigned j = 0; j < DIMENSIONS; ++j)
 140                         {
 141                                 token = strtok(NULL, ",; ");
 142                                 a->x[j] = atof(token);
 143                         }
 144                         for (unsigned j = 0; j < DIMENSIONS; ++j)
 145                         {
 146                                 token = strtok(NULL, ",; ");
 147                                 a->v[j] = atof(token);
 148                         }
 149                         //Body_Print(a);
 150                 }
 151         }
 152
 153         //printf("--------------------------------------------------------------\n");
 154
 155         fclose(file);
 156 }
 157
 158 /**
 159  * Cleans up the universe by freeing all memory
 160  */
 161 void Universe_Cleanup()
 162 {
 163         free(universe.body);
 164
 165 }
 166
 167 /**
 168  * Thread - Continuously computes steps for a system of bodies
 169  */
 170 void * Compute_Thread(void * arg)
 171 {
 172
 173         System * s = (System*)(arg); //cast argument to a System*
 174         // First set up the "sub_system" array
 175         #if NUM_THREADS > 1
 176                 unsigned bodies_per_system = (s->N) / NUM_THREADS;
 177                 unsigned remainder = (s->N) % NUM_THREADS;
 178                 for (unsigned i = 0; i < NUM_THREADS; ++i)
 179                 {
 180                         sub_system[i].body = (s->body)+(i*bodies_per_system);
 181                         sub_system[i].N = bodies_per_system;
 182                         if (i == NUM_THREADS - 1)
 183                                 sub_system[i].N += remainder;
 184                 }
 185
 186                 pthread_attr_t attr; //thread attribute for the workers
 187                 pthread_attr_init(&attr);
 188                 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_DETACHED);
 189
 190         #endif
 191         while (true)
 192         {
 193
 194                 if (runstate != RUN) pthread_exit(NULL); //Check whether the main thread wants to exit
 195
 196
 197
 198                 #if NUM_THREADS <= 1
 199                         //No worker threads; do everything in this thread
 200                         System_Forces(s, s);
 201                         System_Positions(s);
 202                 #else
 203
 204
 205
 206
 207                 worker_count = NUM_THREADS; //All threads working
 208
 209                 //Compute forces
 210                 for (unsigned i = 0; i < NUM_THREADS; ++i)
 211                 {
 212                         if (pthread_create(worker_thread+i, &attr, Force_Thread, (void*)(sub_system+i)) != 0)
 213                         {
 214                                 perror("In compute thread, couldn't create worker thread (force)");
 215                                 QuitProgram(true);
 216                                 pthread_exit(NULL);
 217                         }
 218                 }
 219
 220
 221
 222                 //Wait for forces to be computed
 223                 pthread_mutex_lock(&mutex_workers);
 224                 while (worker_count > 0)
 225                         pthread_cond_wait(&workers_done_cv, &mutex_workers);
 226                 pthread_mutex_unlock(&mutex_workers);
 227
 228                 worker_count = NUM_THREADS; //All threads working
 229
 230                 //Compute positions
 231                 for (unsigned i = 0; i < NUM_THREADS; ++i)
 232                 {
 233                         if (pthread_create(worker_thread+i, &attr, Position_Thread, (void*)(sub_system+i)) != 0)
 234                         {
 235                                 perror("In compute thread, couldn't create worker thread (position)");
 236                                 QuitProgram(true);
 237                                 pthread_exit(NULL);
 238                         }
 239                 }
 240
 241                 //Wait for positions to be computed
 242                 pthread_mutex_lock(&mutex_workers);
 243                 while (worker_count > 0)
 244                         pthread_cond_wait(&workers_done_cv, &mutex_workers);
 245                 pthread_mutex_unlock(&mutex_workers);
 246
 247
 248                 #endif
 249
 250         }
 251 }
 252
 253 void * Force_Thread(void * s)
 254 {
 255
 256         System_Forces((System*)s, &universe);
 257         pthread_mutex_lock(&mutex_workers);
 258         worker_count -= 1;
 259         if (worker_count == 0)
 260         {
 261                 pthread_cond_signal(&workers_done_cv);
 262         }
 263         pthread_mutex_unlock(&mutex_workers);
 264         return NULL;
 265 }
 266
 267 void * Position_Thread(void * s)
 268 {
 269
 270         System_Positions((System*)s);
 271         pthread_mutex_lock(&mutex_workers);
 272         worker_count -= 1;
 273         if (worker_count == 0)
 274         {
 275                 pthread_cond_signal(&workers_done_cv);
 276         }
 277         pthread_mutex_unlock(&mutex_workers);
 278         return NULL;
 279 }
 280
 281 /**
 282  * Child threads can call this to signal the main thread to quit the program
 283  * The main thread also uses this to tell child threads that the program is quitting
 284  *      Note that the function doesn't call exit(), that is done by the main thread
 285  */
 286 void QuitProgram(bool error)
 287 {
 288
 289         pthread_mutex_lock(&mutex_runstate);
 290         if (error)
 291                 runstate = QUIT_ERROR;
 292         else
 293                 runstate = QUIT;
 294         pthread_mutex_unlock(&mutex_runstate);
 295 }
 296
 297