3 * @purpose Barnes Hut OpenMP version of N-Body simulator, implementation
4 * @author Sam Moore (20503628) - 2012
8 #include "../single-thread/nbody.c" // Include original code
9 #include "graphics.h" // For Graphics_Run function only
15 omp_lock_t runstate_lock;
18 omp_lock_t nested_lock;
24 * @function Simulation_Run
25 * @purpose Start the simulation
26 * @param argc, argv - Passed to Graphics_Run if needed, otherwise unused
28 void Simulation_Run(int argc, char ** argv)
31 omp_init_lock(&runstate_lock);
33 omp_init_lock(&nested_lock);
37 if (!options.pedantic_graphics)
39 fprintf(stderr, "Warning: Using pedantic graphics to ensure alternate reality matches our own (as closely as possible).\n");
40 options.pedantic_graphics = true;
42 #endif //DUAL_UNIVERSE
45 if (omp_get_nested() == 0)
47 fprintf(stderr, "Couldn't enable nested parallelism\n.");
51 if (options.num_threads == 0)
53 options.num_threads = omp_get_max_threads();
54 if (options.draw_graphics)
55 options.num_threads -= 1;
57 if (options.nested_threads == 0)
59 options.nested_threads = omp_get_max_threads() - options.num_threads;
60 if (options.nested_threads == 0)
61 options.nested_threads = 1;
65 System_Init(&alternate_universe, options.input);
66 #endif //DUAL_UNIVERSE
68 if (options.draw_graphics)
70 Graphics_Run(argc, argv);
71 if (options.pedantic_graphics == false)
76 #pragma omp parallel num_threads(2)
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
83 //printf("Thread %d gets graphics\n", omp_get_thread_num());
85 QuitProgram(false); // Program gets to here if the window is quit
87 else // #pragma omp section
89 //printf("Thread %d gets computations\n", omp_get_thread_num());
99 omp_destroy_lock(&runstate_lock);
101 omp_destroy_lock(&nested_lock);
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...
118 if (omp_get_nested() == 0)
120 fprintf(stderr, "Couldn't enable nested parallelism\n.");
124 omp_set_num_threads(options.num_threads);
125 if (omp_get_max_threads() != options.num_threads)
127 fprintf(stderr, "Couldn't use %d threads!\n", options.num_threads);
131 parent = Node_Create(NULL);
133 Generate_Tree(&universe, parent);
140 while (!ExitCondition())
144 Node_CalculateMass(parent);
147 #pragma omp for // The easiest way to parallelise it
148 for (unsigned a = 0; a < universe.N; ++a)
150 for (unsigned i = 0; i < DIMENSIONS; ++i)
151 universe.body[a].F[i] = 0;
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);
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)
163 Body_Position(universe.body+a);
166 #pragma omp single // Ideally the first thread out should do this
168 Update_Tree(&universe, parent);
171 // I haven't noticed any difference in performance.
172 // I'm just showing off.
175 // Only the master may do stuff with graphics
178 if (options.draw_graphics && options.pedantic_graphics)
185 // Someone should do this stuff too
189 if (options.verbosity != 0 && universe.steps % options.verbosity == 0)
195 Node_Destroy(parent);
202 void QuitProgram(bool error)
204 if (runstate == QUIT || runstate == QUIT_ERROR)
206 omp_set_lock(&runstate_lock);
207 runstate = (error) ? QUIT_ERROR : QUIT;
208 omp_unset_lock(&runstate_lock);
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);
218 for (int i = 0; i < alternate_universe.N; ++i)
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
230 // Draw a line from the "alternate" body to the real one
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]);
237 #endif //DUAL_UNIVERSE
242 /* // Uncomment to draw cubes
243 glPolygonMode( GL_FRONT_AND_BACK, GL_LINE );
246 glPolygonMode( GL_FRONT_AND_BACK, GL_FILL);