--- /dev/null
+/**
+ * @file nbody.c
+ * @purpose Barnes Hut OpenMP version of N-Body simulator, implementation
+ * @author Sam Moore (20503628) - 2012
+ */
+
+#include "nbody.h"
+#include "../single-thread/nbody.c" // Include original code
+#include "graphics.h" // For Graphics functions
+
+#include "tree.h"
+
+
+
+omp_lock_t runstate_lock;
+
+#ifdef USE_THREADS
+omp_lock_t nested_lock;
+#endif //USE_THREADS
+
+Node * parent = NULL;
+
+/**
+ * @function Simulation_Run
+ * @purpose Start the simulation
+ * @param argc, argv - Passed to Graphics_Run if needed, otherwise unused
+ */
+void Simulation_Run(int argc, char ** argv)
+{
+
+ omp_init_lock(&runstate_lock);
+ #ifdef USE_THREADS
+ omp_init_lock(&nested_lock);
+ #endif //USE_THREADS
+
+ #ifdef DUAL_UNIVERSE
+ if (!options.pedantic_graphics)
+ {
+ fprintf(stderr, "Warning: Using pedantic graphics to ensure alternate reality matches our own (as closely as possible).\n");
+ options.pedantic_graphics = true;
+ }
+ #endif //DUAL_UNIVERSE
+
+ omp_set_nested(1);
+ if (omp_get_nested() == 0)
+ {
+ fprintf(stderr, "Couldn't enable nested parallelism\n.");
+ exit(EXIT_FAILURE);
+ }
+
+ if (options.num_threads == 0)
+ {
+ options.num_threads = omp_get_max_threads();
+ if (options.draw_graphics)
+ options.num_threads -= 1;
+ }
+ if (options.nested_threads == 0)
+ {
+ options.nested_threads = omp_get_max_threads() - options.num_threads;
+ if (options.nested_threads == 0)
+ options.nested_threads = 1;
+ }
+
+ #ifdef DUAL_UNIVERSE
+ System_Init(&alternate_universe, options.input);
+ #endif //DUAL_UNIVERSE
+
+ if (options.draw_graphics)
+ {
+ Graphics_Run(argc, argv);
+ if (options.pedantic_graphics == false)
+ {
+
+
+
+ #pragma omp parallel num_threads(2)
+ {
+ // This can't be done with sections!
+ // Because glut is useless, and can only be dealt with in the main thread
+ // I just hope that thread "0" is always the main thread!
+ if (omp_get_thread_num() == 0) // #pragma omp section
+ {
+ //printf("Thread %d gets graphics\n", omp_get_thread_num());
+ glutMainLoop();
+ QuitProgram(false); // Program gets to here if the window is quit
+ }
+ else // #pragma omp section
+ {
+ //printf("Thread %d gets computations\n", omp_get_thread_num());
+ Compute();
+ }
+
+ }
+ return;
+ }
+ }
+ Compute();
+
+ omp_destroy_lock(&runstate_lock);
+ #ifdef USE_THREADS
+ omp_destroy_lock(&nested_lock);
+ #endif //USE_THREADS
+
+
+
+}
+
+/**
+ * @function Compute
+ * @purpose The main computation loop
+ * Uses threading, but only on the first level
+ * I tried to introduce threading to some of the recursive functions in the tree...
+ * But failed
+ */
+void Compute(void)
+{
+ omp_set_nested(1);
+ if (omp_get_nested() == 0)
+ {
+ fprintf(stderr, "Couldn't enable nested parallelism\n.");
+ exit(EXIT_FAILURE);
+ }
+
+ omp_set_num_threads(options.num_threads);
+ if (omp_get_max_threads() != options.num_threads)
+ {
+ fprintf(stderr, "Couldn't use %d threads!\n", options.num_threads);
+ exit(EXIT_FAILURE);
+ }
+
+ parent = Node_Create(NULL);
+
+ Generate_Tree(&universe, parent);
+
+
+
+
+#pragma omp parallel
+{
+ while (!ExitCondition())
+ {
+ #pragma omp single
+ {
+ Node_CalculateMass(parent);
+ }
+
+ #pragma omp for // The easiest way to parallelise it
+ for (unsigned a = 0; a < universe.N; ++a)
+ {
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ universe.body[a].F[i] = 0;
+
+ Node_Forces(universe.body+a, parent, options.theta);
+ //printf("Forces on %p done\n", (void*)(universe.body+a));
+ Body_Velocity(universe.body+a);
+ }
+
+ // put a nowait here so that the first thread out can start updating the tree
+ #pragma omp for nowait
+ for (unsigned a = 0; a < universe.N; ++a)
+ {
+
+ Body_Position(universe.body+a);
+ }
+
+ #pragma omp single // Ideally the first thread out should do this
+ {
+ Update_Tree(&universe, parent);
+ }
+
+ // I haven't noticed any difference in performance.
+ // I'm just showing off.
+ #pragma omp barrier
+
+ // Only the master may do stuff with graphics
+ #pragma omp master
+ {
+ if (options.draw_graphics && options.pedantic_graphics)
+ {
+ glutMainLoopEvent();
+ Graphics_Display();
+ }
+ }
+
+ // Someone should do this stuff too
+ #pragma omp single
+ {
+ universe.steps += 1;
+ if (options.verbosity != 0 && universe.steps % options.verbosity == 0)
+ DisplayStatistics();
+ }
+
+ }
+}
+ Node_Destroy(parent);
+ free(parent);
+
+ QuitProgram(false);
+}
+
+
+void QuitProgram(bool error)
+{
+ if (runstate == QUIT || runstate == QUIT_ERROR)
+ return;
+ omp_set_lock(&runstate_lock);
+ runstate = (error) ? QUIT_ERROR : QUIT;
+ omp_unset_lock(&runstate_lock);
+}
+
+#ifdef DUAL_UNIVERSE
+void BeforeDraw()
+{
+ // Before drawing the main universe... compute and draw the alternate universe
+ // Not threaded, because I am too lazy
+ System_Compute(&alternate_universe);
+
+ for (int i = 0; i < alternate_universe.N; ++i)
+ {
+
+ Body * b = alternate_universe.body+i;
+ Body * original = universe.body+i;
+ //glColor3f(0.0f, b->mass/1e11*100, 0.0f);
+ glColor3f(1.0f, 0.0f, 1.0f);
+ glPushMatrix(); // to save the current matrix
+ glTranslated(scale*b->x[0], scale*b->x[1], scale*b->x[2]);
+ glutSolidSphere (BALL_SIZE, 10, 10);
+ glPopMatrix(); // restore the previous matrix
+
+ // Draw a line from the "alternate" body to the real one
+ glBegin(GL_LINES);
+ glVertex3f(scale*b->x[0], scale*b->x[1], scale*b->x[2]);
+ glVertex3f(scale*original->x[0], scale*original->x[1], scale*original->x[2]);
+ glEnd();
+ }
+}
+#endif //DUAL_UNIVERSE
+
+#ifndef AfterDraw
+void AfterDraw()
+{
+ /* // Uncomment to draw cubes
+ glPolygonMode( GL_FRONT_AND_BACK, GL_LINE );
+ Draw_Node(parent);
+ glFlush();
+ glPolygonMode( GL_FRONT_AND_BACK, GL_FILL);
+ */
+}
+#endif //AfterDraw
git.ucc.asn.au / matches / honours.git / blobdiff
UCC git Repository :: git.ucc.asn.au