course/semester2/pprog/assignment1/nbody-bh/nbody.c

   1 /**
   2  * @file nbody.c
   3  * @purpose Barnes Hut OpenMP version of N-Body simulator, implementation
   4  * @author Sam Moore (20503628) - 2012
   5  */
   6
   7 #include "nbody.h"
   8 #include "../single-thread/nbody.c" // Include original code
   9 #include "graphics.h" // For Graphics functions
  10
  11 #include "tree.h"
  12
  13
  14
  15 omp_lock_t runstate_lock;
  16
  17 #ifdef USE_THREADS
  18 omp_lock_t nested_lock;
  19 #endif //USE_THREADS
  20
  21 Node * parent = NULL;
  22
  23 /**
  24  * @function Simulation_Run
  25  * @purpose Start the simulation
  26  * @param argc, argv - Passed to Graphics_Run if needed, otherwise unused
  27  */
  28 void Simulation_Run(int argc, char ** argv)
  29 {
  30
  31         omp_init_lock(&runstate_lock);
  32         #ifdef USE_THREADS
  33         omp_init_lock(&nested_lock);
  34         #endif //USE_THREADS
  35
  36         #ifdef DUAL_UNIVERSE
  37         if (!options.pedantic_graphics)
  38         {
  39                 fprintf(stderr, "Warning: Using pedantic graphics to ensure alternate reality matches our own (as closely as possible).\n");
  40                 options.pedantic_graphics = true;
  41         }
  42         #endif //DUAL_UNIVERSE
  43
  44         omp_set_nested(1);
  45         if (omp_get_nested() == 0)
  46         {
  47                 fprintf(stderr, "Couldn't enable nested parallelism\n.");
  48                 exit(EXIT_FAILURE);
  49         }
  50
  51         if (options.num_threads == 0)
  52         {
  53                 options.num_threads = omp_get_max_threads();
  54                 if (options.draw_graphics)
  55                         options.num_threads -= 1;
  56         }
  57         if (options.nested_threads == 0)
  58         {
  59                 options.nested_threads = omp_get_max_threads() - options.num_threads;
  60                 if (options.nested_threads == 0)
  61                         options.nested_threads = 1;
  62         }
  63
  64         #ifdef DUAL_UNIVERSE
  65                 System_Init(&alternate_universe, options.input);
  66         #endif //DUAL_UNIVERSE
  67
  68         if (options.draw_graphics)
  69         {
  70                 Graphics_Run(argc, argv);
  71                 if (options.pedantic_graphics == false)
  72                 {
  73
  74
  75
  76                         #pragma omp parallel num_threads(2)
  77                         {
  78                                 // This can't be done with sections!
  79                                 // Because glut is useless, and can only be dealt with in the main thread
  80                                 // I just hope that thread "0" is always the main thread!
  81                                 if (omp_get_thread_num() == 0) // #pragma omp section
  82                                 {
  83                                         //printf("Thread %d gets graphics\n", omp_get_thread_num());
  84                                         glutMainLoop();
  85                                         QuitProgram(false); // Program gets to here if the window is quit
  86                                 }
  87                                 else // #pragma omp section
  88                                 {
  89                                         //printf("Thread %d gets computations\n", omp_get_thread_num());
  90                                         Compute();
  91                                 }
  92
  93                         }
  94                         return;
  95                 }
  96         }
  97         Compute();
  98
  99         omp_destroy_lock(&runstate_lock);
 100         #ifdef USE_THREADS
 101         omp_destroy_lock(&nested_lock);
 102         #endif //USE_THREADS
 103
 104
 105
 106 }
 107
 108 /**
 109  * @function Compute
 110  * @purpose The main computation loop
 111  *      Uses threading, but only on the first level
 112  *      I tried to introduce threading to some of the recursive functions in the tree...
 113  *      But failed
 114  */
 115 void Compute(void)
 116 {
 117         omp_set_nested(1);
 118         if (omp_get_nested() == 0)
 119         {
 120                 fprintf(stderr, "Couldn't enable nested parallelism\n.");
 121                 exit(EXIT_FAILURE);
 122         }
 123
 124         omp_set_num_threads(options.num_threads);
 125         if (omp_get_max_threads() != options.num_threads)
 126         {
 127                 fprintf(stderr, "Couldn't use %d threads!\n", options.num_threads);
 128                 exit(EXIT_FAILURE);
 129         }
 130
 131         parent = Node_Create(NULL);
 132
 133         Generate_Tree(&universe, parent);
 134
 135
 136
 137
 138 #pragma omp parallel
 139 {
 140         while (!ExitCondition())
 141         {
 142                 #pragma omp single
 143                 {
 144                         Node_CalculateMass(parent);
 145                 }
 146
 147                 #pragma omp for // The easiest way to parallelise it
 148                 for (unsigned a = 0; a < universe.N; ++a)
 149                 {
 150                         for (unsigned i = 0; i < DIMENSIONS; ++i)
 151                                 universe.body[a].F[i] = 0;
 152
 153                         Node_Forces(universe.body+a, parent, options.theta);
 154                         //printf("Forces on %p done\n", (void*)(universe.body+a));
 155                         Body_Velocity(universe.body+a);
 156                 }
 157
 158                 // put a nowait here so that the first thread out can start updating the tree
 159                 #pragma omp for nowait
 160                 for (unsigned a = 0; a < universe.N; ++a)
 161                 {
 162
 163                         Body_Position(universe.body+a);
 164                 }
 165
 166                 #pragma omp single // Ideally the first thread out should do this
 167                 {
 168                         Update_Tree(&universe, parent);
 169                 }
 170
 171                 // I haven't noticed any difference in performance.
 172                 // I'm just showing off.
 173                 #pragma omp barrier
 174
 175                 // Only the master may do stuff with graphics
 176                 #pragma omp master
 177                 {
 178                         if (options.draw_graphics && options.pedantic_graphics)
 179                         {
 180                                 glutMainLoopEvent();
 181                                 Graphics_Display();
 182                         }
 183                 }
 184
 185                 // Someone should do this stuff too
 186                 #pragma omp single
 187                 {
 188                         universe.steps += 1;
 189                         if (options.verbosity != 0 && universe.steps % options.verbosity == 0)
 190                                 DisplayStatistics();
 191                 }
 192
 193         }
 194 }
 195         Node_Destroy(parent);
 196         free(parent);
 197
 198         QuitProgram(false);
 199 }
 200
 201
 202 void QuitProgram(bool error)
 203 {
 204         if (runstate == QUIT || runstate == QUIT_ERROR)
 205                 return;
 206         omp_set_lock(&runstate_lock);
 207         runstate = (error) ? QUIT_ERROR : QUIT;
 208         omp_unset_lock(&runstate_lock);
 209 }
 210
 211 #ifdef DUAL_UNIVERSE
 212 void BeforeDraw()
 213 {
 214         // Before drawing the main universe... compute and draw the alternate universe
 215         // Not threaded, because I am too lazy
 216         System_Compute(&alternate_universe);
 217
 218         for (int i = 0; i < alternate_universe.N; ++i)
 219         {
 220
 221                 Body * b = alternate_universe.body+i;
 222                 Body * original = universe.body+i;
 223                 //glColor3f(0.0f, b->mass/1e11*100, 0.0f);
 224                 glColor3f(1.0f, 0.0f, 1.0f);
 225                 glPushMatrix(); // to save the current matrix
 226                 glTranslated(scale*b->x[0], scale*b->x[1], scale*b->x[2]);
 227                 glutSolidSphere (BALL_SIZE, 10, 10);
 228                 glPopMatrix(); // restore the previous matrix
 229
 230                 // Draw a line from the "alternate" body to the real one
 231                 glBegin(GL_LINES);
 232                 glVertex3f(scale*b->x[0], scale*b->x[1], scale*b->x[2]);
 233                 glVertex3f(scale*original->x[0], scale*original->x[1], scale*original->x[2]);
 234                 glEnd();
 235         }
 236 }
 237 #endif //DUAL_UNIVERSE
 238
 239 #ifndef AfterDraw
 240 void AfterDraw()
 241 {
 242         /* // Uncomment to draw cubes
 243         glPolygonMode( GL_FRONT_AND_BACK, GL_LINE );
 244         Draw_Node(parent);
 245         glFlush();
 246         glPolygonMode( GL_FRONT_AND_BACK, GL_FILL);
 247         */
 248 }
 249 #endif //AfterDraw