-/**
- * @function Split_System
- * @purpose Helper to divide one system into an array of systems
- * Each sub system will have N = (s->N / n) bodies in it
- * @param s - The original system (typically &universe)
- * @param n - The number of sub systems in the array
- *
- * WARNING: It is the caller's responsibility to free() the returned array
- */
-System * Split_System(System * s, unsigned n)
-{
- System * result = (System*)(calloc(n, sizeof(System)));
- if (result == NULL)
- {
- perror("Couldn't create array of sub systems");
- QuitProgram(true);
- pthread_exit(NULL);
- }
- unsigned n_per_system = (s->N) / n;
- unsigned remainder = (s->N) % n;
-
- for (unsigned i = 0; i < n; ++i)
- {
- result[i].N = n_per_system;
- if (i == n-1)
- result[i].N += remainder;
- result[i].body = (s->body) + (n_per_system * i);
- result[i].steps = 0;
- }
- return result;
-}
/**
* @function Allocate_Threads
System * s = (System*)(arg);
s->steps += 1; //Increment number of steps computed
- if (options.verbosity != 0 && s->steps % options.verbosity == 1)
+ if (options.verbosity != 0 && s->steps % options.verbosity == 0)
DisplayStatistics();
--- /dev/null
+/**
+ * @file nbody.c
+ * @purpose Implentation of multi-threaded N-Body simulator using pthreads
+ * @author Sam Moore (20503628) - 2012
+ */
+
+#include "nbody.h" // Declarations
+#include "../single-thread/nbody.c" // Include all functions from the single threaded version
+#include <assert.h>
+#include "graphics.h" // For declaration of Graphics_Run only
+#include "barrier.c"
+
+// --- Variable declarations --- //
+
+pthread_t compute_thread; // The thread responsible for computations; it spawns worker threads (* terms and conditions apply)
+
+pthread_t * worker_thread = NULL; //Array of worker threads responsible for Force and Position updates
+System * sub_system = NULL; //Array of Systems used to divide up the main "universe" System for worker threads
+System * nested_sub_system = NULL; //Array of Systems for division of "universe" for the nested worker threads
+
+pthread_mutex_t mutex_runstate; // Mutex around the runstate
+
+
+
+
+Barrier force_barrier; // I laughed at this variable name. A bit sad really.
+Barrier position_barrier;
+
+
+Barrier graphics_barrier;
+
+pthread_mutex_t mutex_threads_running;
+int threads_running = 0;
+
+/**
+ * @function Compute_Thread
+ * @purpose Thread - Continuously computes steps for a system of bodies. Seperate from graphics functions.
+ * Spawns worker threads to divide up computation.
+ * @param arg - Can be cast to the System for which computations are performed
+ *
+ * NOTE:
+ * This will always be (void*)(&universe) where universe is the global System variable that has every Body in it.
+ * But I don't like global variables. And since the argument had to be passed, I thought I might as well use it.
+ * That way, when I change "universe" to "solar_system", I only have to change the argument where this is called, not all through it.
+ * Find and replace? Who uses that!?
+ */
+
+
+void * Compute_Thread(void * arg)
+{
+
+ System * s = (System*)(arg); //cast argument to a System*
+
+
+
+ // If no number of threads provided, use the default value, unless someone (me) changed that to a stupid value
+ if (options.num_threads <= 0)
+ options.num_threads = (DEFAULT_WORKING_THREADS > 1) ? DEFAULT_WORKING_THREADS : 1; // Fear the ternary operator!
+
+ if (options.nested_threads <= 0)
+ options.nested_threads = 1;
+
+ // Do a sanity check; there is no point spawning more threads than bodies.
+ if (options.num_threads > s->N)
+ {
+ fprintf(stderr,
+ "Warning: Using %u threads instead of %u specified, because there are only %u bodies to simulate!\n",
+ s->N, options.num_threads, s->N);
+ options.num_threads = s->N;
+ }
+
+ // Initialise the barriers ("shields up!")
+ Barrier_Init(&force_barrier, options.num_threads);
+ Barrier_Init(&position_barrier, options.num_threads);
+ Barrier_Init(&graphics_barrier, 1);
+
+
+
+ if (options.num_threads > 1) // If we require additional worker threads...
+ {
+ // Allocate worker threads and sub systems
+
+ sub_system = Split_System(&universe, options.num_threads);
+
+ if (options.nested_threads > 1)
+ nested_sub_system = Split_System(&universe, options.nested_threads);
+
+ #ifdef PERSISTENT_THREADS // Code for the smart way of doing it (spawn threads once, keep them running)
+ worker_thread = Allocate_Threads(options.num_threads-1);
+ // Spawn a bunch of worker threads, and let them all do their thing.
+ // Note the "-1". Because this thread can do work too!
+ for (unsigned i = 0; i < options.num_threads-1; ++i)
+ {
+ if (pthread_create(worker_thread+i, NULL, Worker_Thread, (void*)(sub_system+i)) != 0)
+ {
+ perror("In compute thread, couldn't create worker thread");
+ QuitProgram(true);
+ pthread_exit(NULL);
+ }
+ }
+ Worker_Thread((void*)(sub_system+options.num_threads-1)); // This thread becomes a worker thread
+ #else
+ worker_thread = Allocate_Threads(options.num_threads);
+ #endif //PERSISTENT_THREADS
+
+ }
+ #ifdef PERSISTENT_THREADS
+ else // We only require one worker thread...
+ {
+ // So just do all computations in this thread
+ while (!ExitCondition())
+ {
+
+ System_Forces(s, s);
+ // If required, wait for graphics to finish drawing
+ if (options.draw_graphics && options.pedantic_graphics)
+ Barrier_Wait(&graphics_barrier);
+ System_Positions(s);
+ StepFunction(s);
+
+ if (options.draw_graphics && options.pedantic_graphics)
+ Barrier_ForceExit(&position_barrier); //Make the graphics continue
+ }
+ QuitProgram(false);
+ pthread_exit(NULL);
+ }
+
+ #else // Code for the stupid way of doing it (respawn threads each step)
+ // (ie: The way I immediately implemented and didn't realise was stupid until someone told me)
+
+
+
+
+ // Run until we can't run anymore
+ while (!ExitCondition())
+ {
+
+
+ if (options.num_threads <= 1) // If there is only 1 worker thread...
+ {
+ // Just do everything in this thread
+ System_Forces(s, s);
+ // If required, wait for graphics to finish drawing
+ if (options.draw_graphics && options.pedantic_graphics)
+ Barrier_Join(&graphics_barrier);
+
+ System_Positions(s);
+ StepFunction(s);
+
+ if (options.draw_graphics && options.pedantic_graphics)
+ Barrier_Join(&position_barrier); //Make the graphics continue
+ continue;
+ }
+
+
+
+ //Compute forces by spawning threads, each thread gets a sub system
+ for (unsigned i = 0; i < options.num_threads; ++i)
+ {
+ if (pthread_create(worker_thread+i, NULL, Force_Thread, (void*)(sub_system+i)) != 0)
+ {
+ perror("In compute thread, couldn't create worker thread (force)");
+ QuitProgram(true);
+ pthread_exit(NULL);
+ }
+
+ }
+
+ for (unsigned i = 0; i < options.num_threads; ++i)
+ pthread_join(worker_thread[i], NULL);
+
+
+ // If required, wait for graphics to finish drawing
+ if (options.draw_graphics && options.pedantic_graphics)
+ Barrier_Wait(&graphics_barrier);
+
+
+ Barrier_Enter(&position_barrier);
+
+ //Compute positions by spawning a bunch of threads to do it
+ for (unsigned i = 0; i < options.num_threads; ++i)
+ {
+ if (pthread_create(worker_thread+i, NULL, Position_Thread, (void*)(sub_system+i)) != 0)
+ {
+ perror("In compute thread, couldn't create worker thread (position)");
+ QuitProgram(true);
+ pthread_exit(NULL);
+ }
+ }
+
+
+ for (unsigned i = 0; i < options.num_threads; ++i)
+ pthread_join(worker_thread[i], NULL);
+
+
+ StepFunction(s); // Execute single threaded stuff
+
+
+ }
+ QuitProgram(false);
+ #endif //PERSISTENT_THREADS
+ return NULL;
+}
+
+/**
+ * @function BeforeDraw
+ * @purpose Called in graphics thread before the draw loop
+ * When --pedantic-graphics enabled, will wait for position computations to finish before drawing
+ * Otherwise does nothing
+ *
+ * This originally seemed like a good place to put the code now in StepFunction(), since only one thread runs this
+ * But then I realised that the graphics might be disabled,
+ * and there was no point having a thread that only existed to call that code.
+ *
+ * So I changed it to the horrible solution that I currently have.
+ */
+void BeforeDraw()
+{
+
+ //printf("BEFORE DRAW\n");
+ if (!options.pedantic_graphics)
+ return;
+
+ //printf("Graphics thread waits on position barrier\n");
+ Barrier_Wait(&position_barrier);
+ //printf("\tGraphics thread wakes up\n");
+ Barrier_Enter(&graphics_barrier);
+}
+
+/**
+ * @function AfterDraw
+ * @purpose Called in graphics thread after the draw loop
+ * When --pedantic-graphics is supplied, will signal computation threads that drawing is finished
+ * So that positions can be safely altered
+ * Otherwise does nothing
+ */
+void AfterDraw()
+{
+ //universe.steps += 1;
+ if (!options.pedantic_graphics)
+ return;
+ Barrier_Join(&graphics_barrier);
+
+}
+
+#ifdef PERSISTENT_THREADS
+
+/**
+ * @function Worker_Thread
+ * @purpose Thread - A self contained worker thread to compute a particular sub system of bodies
+ *
+ * This is the "smart" way to do it, because threads are only created once, and compute both force and position.
+ * The greatest difficulty with pthreads is getting a *single* thread from the team to execute certain code
+ * (ie: The stuff in StepFunction()).
+ * With the "continuously respawning threads of stupidity" approach,
+ * because there is one "master" thread (not necessarilly the main thread... don't get confused now)
+ * to keep respawning the workers, the single threaded code can just be executed in the master thread.
+ *
+ * With this approach, I have created a hacky solution so that the *last* thread to leave the position barrier gets to call StepFunction.
+ *
+ */
+void * Worker_Thread(void * arg)
+{
+ System * s = (System*)(arg); // This is mainly to save typing the RHS a lot of times
+
+ // Each thread runs until the whole program is supposed to end
+ while (!ExitCondition())
+ {
+
+
+ System_Forces(s, &universe); // Each thread computes the forces for its share of bodies
+
+ // Do not confuse with "Barrier_Wait".
+ // Barrier_Wait does not affect the barrier; it just waits for it
+ // Barrier_Join actively updates the state of the barrier, and wakes up sleeping threads if required.
+
+ Barrier_Join(&force_barrier); // All threads must reach here before moving on.
+ if (ExitCondition()) return NULL;
+
+
+ //fprintf(stderr,"Thread %p - force barrier finished\n", arg);
+ //printf("Computed ALL forces\n");
+
+
+ // If required, wait for the graphics to finish drawing stuff
+ if (options.draw_graphics && options.pedantic_graphics)
+ {
+ //printf("Worker %p waits on graphics barrier\n", arg);
+ Barrier_Wait(&graphics_barrier);
+ //printf("\tWorker %p wakes up after graphics barrier\n", arg);
+ if (ExitCondition()) return NULL;
+ }
+
+
+
+ Barrier_Enter(&position_barrier);
+ System_Positions(s); // Each thread updates the positions for its share of bodies
+
+
+ // Barrier_JoinCall behaves in the same way as Barrier_Join, except the *last* thread
+ // (ie: the one that wakes up the others) also calls the function with arguments given.
+ Barrier_JoinCall(&position_barrier, StepFunction, (void*)(&universe));
+ if (ExitCondition()) return NULL;
+ //Barrier_Join(&position_barrier);
+
+ // All threads have computed positions, and *one* thread calls StepFunction()
+
+ }
+ QuitProgram(false); // Set the run state of the program
+ return NULL;
+}
+
+#endif //PERSISTENT_THREADS
+
+/**
+ * @function Force_Thread
+ * @purpose Thread - Calculates the forces on objects in a System
+ * @param s - Can be cast to the System*
+ */
+void * Force_Thread(void * s)
+{
+ //System_ForcePair * pair = (System_ForcePair*)s;
+
+
+ System_Forces(s, &universe); //Simple wrapper
+ //printf("Force_Thread waits\n");
+
+ return NULL;
+}
+
+
+/**
+ * @function Position_Thread
+ * @purpose Thread - Calculates the positions of objects in a System
+ * @param s - Can be cast to the System*
+ */
+void * Position_Thread(void * s)
+{
+
+ System_Positions((System*)s); // Simple wrapper
+ Barrier_Join(&position_barrier); // This needed so that graphics will wait
+
+ return NULL;
+}
+
+/**
+ * @function QuitProgram
+ * @purpose This function can be called in any thread to signal all threads to exit
+ * Repeated calls to this function have no effect
+ *
+ * All threads periodically call ExitCondition(), which will return true if the program should exit.
+ * One (not the only way) to return true is if this function has been called.
+ * Threads will call this function if they detect ExitCondition() is true. Only the first call has any effect.
+ */
+inline void QuitProgram(bool error)
+{
+ //If already quitting, don't do anything
+ if (runstate == QUIT || runstate == QUIT_ERROR)
+ return;
+
+
+
+ // set the runstate (checked in ExitCondition())
+
+ pthread_mutex_lock(&mutex_runstate); // aquire mutex
+ if (error)
+ runstate = QUIT_ERROR; // Program is exiting due to an error
+ else
+ runstate = QUIT; // Program is exiting naturally
+ pthread_mutex_unlock(&mutex_runstate); //release mutex
+
+
+}
+
+/**
+ * @function Thread_Cleanup
+ * @purpose Will be called in the *main* thread when exit() is called
+ * Ensures working threads will exit, and waits for them to finish.
+ * Then waits for them to finish.
+ * Also frees memory associated with the worker threads.
+ */
+void Thread_Cleanup(void)
+{
+
+
+ // Threads recheck the exit condition whenever they leave a barrier.
+ // These calls will stop any threads waiting forever in a barrier for threads that exited before getting to the barrier.
+ Barrier_ForceExit(&force_barrier);
+ Barrier_ForceExit(&position_barrier);
+
+
+ if (options.draw_graphics) // If the graphics are enabled...
+ {
+ // Then there is a computation thread, since graphics are done in the main thread
+ pthread_join(compute_thread, NULL);
+ }
+
+ #ifdef PERSISTENT_THREADS
+ for (unsigned i = 0; i < options.num_threads-1; ++i)
+ {
+ pthread_join(worker_thread[i], NULL);
+ }
+ #else
+ // All other worker threads (if they were spawned) are terminated in Compute_Thread
+ #endif //PERSISTENT_THREADS
+
+ // Scary memory management here.
+ if (worker_thread != NULL)
+ free(worker_thread);
+ if (sub_system != NULL)
+ free(sub_system);
+ worker_thread = NULL;
+ sub_system = NULL;
+
+}
+
+
+/**
+ * @function Simulation_Run
+ * @purpose Initialise and start the simulation. Will be called in the main thread.
+ * Replaces the single-threaded macro that does nothing, and sets up the graphics and computation threads
+ * @param argc - Number of arguments - Passed to Graphics_Run if needed
+ * @param argv - Argument strings - Passed to Graphics_Run if needed
+ */
+void Simulation_Run(int argc, char ** argv)
+{
+ atexit(Thread_Cleanup);
+ if (options.draw_graphics) // The graphics are enabled
+ {
+ // I have chosen to do graphics in the main thread in this case.
+ // A *single* seperate thread is spawned here to do computations.
+ // This computation thread will spawn any additional worker threads required.
+ if (pthread_create(&compute_thread, NULL, Compute_Thread, (void*)&universe) != 0)
+ {
+ perror("Error creating compute thread");
+ exit(EXIT_FAILURE);
+ }
+
+ // This is run in the main thread
+ // It is effectively the graphics initialisation, followed by the glut loop
+ Graphics_Run(argc, argv);
+
+ // The main thread reaches here after leaving the glut loop when ExitCondition() returns true.
+
+ QuitProgram(false);
+
+ exit(EXIT_SUCCESS); // This is the main thread; use exit()
+
+ }
+ else //The graphics are disabled
+ {
+ // If graphics are disabled, there is no point spawning an extra thread.
+ // In this case, the *main* thread starts computations.
+ // Note that it will probably spawn additional worker threads (unless options.num_threads <= 1)
+ Compute_Thread((void*)(&universe));
+ QuitProgram(false);
+ exit(EXIT_SUCCESS);
+ }
+}
+
+
+
+
+
+/**
+ * @function Allocate_Threads
+ * @purpose Helper function to allocate an array of pthread_t objects
+ * Handles all the pointless, er, "important" error checking that should be done
+ * @param n - Number of threads in the array
+ *
+ * WARNING: Remember to free() the array!!!
+ */
+pthread_t * Allocate_Threads(unsigned n)
+{
+ pthread_t * result = (pthread_t*)(calloc(n, sizeof(pthread_t)));
+ if (result == NULL)
+ {
+ perror("Unable to allocate memory for threads");
+ QuitProgram(true);
+ pthread_exit(NULL);
+ }
+ return result;
+}
+
+/**
+ * @function StepFunction
+ * @purpose Helper to perform stuff in a single thread every step, after position computations are done
+ * The reason this has void* all over the place is so that I can pass the function pointer (to horrible dragons and fiendish demons).
+ * @param arg - Can be cast to System* for which steps are to be updated
+ * Will always be (void*)(&universe). But I have been brainwashed into the "global variables are baaaaad" philosophy.
+ * @returns arg
+ */
+void * StepFunction(void * arg)
+{
+ //fprintf(stderr, "StepFunction called\n");
+ System * s = (System*)(arg);
+ s->steps += 1; //Increment number of steps computed
+
+ if (options.verbosity != 0 && s->steps % options.verbosity == 1)
+ DisplayStatistics();
+
+
+ return arg;
+}
+
-System * Split_System(System * s, unsigned n); // Splits one system into a number of other systems, returns an array of size n
+
pthread_t * Allocate_Threads(unsigned n); // Allocates space for threads - handles errors
#ifdef PERSISTENT_THREADS
--- /dev/null
+#ifndef _NBODY_MTHREAD_H
+#define _NBODY_MTHREAD_H
+
+
+#include "../single-thread/nbody.h" //Use original simulation code
+#include <pthread.h>
+#include "barrier.h" // Barriers
+
+#undef SINGLE_THREADED
+#define PTHREADED
+
+#define DEFAULT_WORKING_THREADS 2
+
+//#define PERSISTENT_THREADS //If defined, threads will not be continually destroyed and then respawned
+
+
+
+//Undefine default macros, replace with functions
+#undef Simulation_Run
+void Simulation_Run(int argc, char ** argv);
+#undef QuitProgram
+void QuitProgram(bool error);
+
+#undef BeforeDraw
+void BeforeDraw();
+#undef AfterDraw
+void AfterDraw();
+
+
+void * Compute_Thread(void * system); //Thread - Continuously perform computations for a System of bodies. May spawn additional worker threads.
+
+
+
+
+System * Split_System(System * s, unsigned n); // Splits one system into a number of other systems, returns an array of size n
+pthread_t * Allocate_Threads(unsigned n); // Allocates space for threads - handles errors
+
+#ifdef PERSISTENT_THREADS
+void * Worker_Thread(void * arg);
+#endif //PERSISTENT_THREADS
+void * Force_Thread(void * system); //Thread - Compute forces for all objects in a system
+void * Position_Thread(void * system); //Thread - Compute positions for all objects in a system
+
+void Thread_Cleanup(void); //Called at program exit to safely join computation thread
+
+
+
+void * StepFunction(void * s); //Function called at the end of every step
+
+#endif //_NBODY_MTHREAD_H
--- /dev/null
+/**
+ * @file graphics.c
+ * @author Sam Moore (20503628) 2012
+ * - Adapted from template program provided by UWA
+ * @purpose Definition of graphics related functions
+ * NOTE: This file is identical for both the single-threaded and multi-threaded versions of the program
+ */
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <time.h> // Needed for clock
+#include <math.h> // Maths functions (sin, cos)
+
+#include "graphics.h" //Function declarations
+#include "nbody.h" //The simulation
+
+
+// --- Variable definitions --- //
+double previousTime, eyeTheta, eyePhi, eyeRho;
+float look[3];
+int windowWidth, windowHeight, upY;
+double scale = 1.0;
+
+#ifdef FLYING_CAMERA
+Camera eye;
+#endif //FLYING_CAMERA
+
+/**
+ * @function Graphics_Run
+ * @purpose Setup and start the graphics engine
+ * @param argc - Number of arguments to main function, passed to glutInit
+ * @param argv - Argument strings for main function, passed to glutInit
+ */
+void Graphics_Run(int argc, char ** argv)
+{
+ if (options.draw_graphics == false)
+ {
+ // This message is left here for when I inevitably accidentally call the function
+ fprintf(stderr, "Graphics_Run should not be called if graphics are disabled!\n");
+ exit(EXIT_FAILURE);
+ }
+
+ glutInit(&argc, argv);
+ glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
+ glutInitWindowSize(WIDTH, HEIGHT);
+ glutInitWindowPosition(POSITION_X, POSITION_Y);
+
+ //Set window title based on version of program
+ #ifdef SINGLE_THREADED
+ glutCreateWindow("N-Body : Single-Threaded");
+ #elif defined PTHREADED
+ glutCreateWindow("N-Body : Multi-Threaded (pthread)");
+ #elif defined OMP_THREADED
+ glutCreateWindow("N-Body : Multi-Threaded (OpenMP)");
+ #else
+ glutCreateWindow("N-Body");
+ #endif
+ glutDisplayFunc(Graphics_Display);
+ glutIdleFunc(Graphics_Display);
+ glutKeyboardFunc(Graphics_Keyboard);
+ glutReshapeFunc(Graphics_Reshape);
+
+
+ glClearColor(1.0,1.0,1.0,0.0);
+ glColor3f(0.0f, 0.0f, 0.0f);
+ glPointSize(POINT_SIZE);
+ glMatrixMode(GL_PROJECTION);
+ glLoadIdentity();
+
+ /*init lighting */
+
+ GLfloat mat_specular[] = { 1.0, 1.0, 1.0, 1.0 };
+ GLfloat mat_shininess[] = { 50.0 };
+ GLfloat light_position[] = { 1.0, 1.0, 0.0, 0.0 };
+ glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
+ glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
+ glLightfv(GL_LIGHT0, GL_POSITION, light_position);
+
+ glColorMaterial(GL_FRONT,GL_DIFFUSE); // Set Color Capability
+
+ glEnable(GL_LIGHTING);
+ glEnable(GL_LIGHT0);
+ glEnable(GL_DEPTH_TEST);
+
+ glEnable(GL_COLOR_MATERIAL); // Enable color
+
+
+ windowWidth = WIDTH;
+ windowHeight = HEIGHT;
+ previousTime = clock();
+ eyeTheta = 0;
+ eyePhi = 0.5 * M_PI;
+ eyeRho = RHO;
+ upY = 1;
+ look[0] = 0;
+ look[1] = 0;
+ look[2] = 0;
+ gluPerspective(VIEW_ANGLE, (double)WIDTH/(double)HEIGHT, WORLD_NEAR, WORLD_FAR);
+
+ #ifdef FLYING_CAMERA
+ eye.x[0] = 0; eye.x[1] = 0; eye.x[2] = 1;
+ eye.y[0] = 0; eye.y[1] = 1; eye.y[2] = 0;
+ eye.z[0] = 1; eye.z[1] = 0; eye.z[2] = 0;
+ eye.p[0] = 0; eye.p[1] = 0; eye.p[2] = -50;
+ #endif //FLYING_CAMERA
+
+
+ //This option must be set, or when glut leaves the main loop. the exit(3) function is called... annoying
+ glutSetOption(GLUT_ACTION_ON_WINDOW_CLOSE, GLUT_ACTION_CONTINUE_EXECUTION);
+
+
+ #ifndef OMP_THREADED
+ glutMainLoop();
+ #endif //OMP_THREADED
+}
+
+/**
+ * @function Graphics_Display
+ * @purpose This function redraws the screen after the positions of particles
+ * have been updated.
+ * It also calls System_Compute, and checks for exit conditions, in the single-threaded version only
+ */
+void Graphics_Display()
+{
+
+ if (options.draw_graphics == false)
+ {
+ // This message is left here for when I inevitably accidentally call the function
+ fprintf(stderr, "Graphics_Display should not be called if graphics are disabled!\n");
+ exit(EXIT_FAILURE);
+ }
+
+ //Check whether the graphics thread should exit for any reason
+ if (ExitCondition())
+ {
+ glutLeaveMainLoop(); // original glut does not have this, which makes "nicely" exiting a bit annoying
+ return;
+ }
+
+ glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+ glMatrixMode(GL_MODELVIEW);
+ glLoadIdentity();
+
+ #ifdef FLYING_CAMERA
+
+ gluLookAt(eye.p[0], eye.p[1], eye.p[2],
+ eye.p[0] + eye.x[0], eye.p[1] + eye.x[1], eye.p[2] + eye.x[2],
+ eye.z[0], eye.z[1], eye.z[2]);
+ #else
+ gluLookAt(eyeRho * sin(eyePhi) * sin(eyeTheta), eyeRho * cos(eyePhi),
+ eyeRho * sin(eyePhi) * cos(eyeTheta),look[0], look[1], look[2], 0, upY, 0);
+ #endif //FLYING_CAMERA
+
+ BeforeDraw(); // Stuff to do before graphics is allowed to draw
+ // Single-thread - perform a computation step, obviously! It's not done anywhere else
+ // Multiple threads - ensure that all body positions are updated (ie: not halfway through step).
+ // (We don't care about the force and velocity variables though)
+
+ for (int i = 0; i < universe.N; ++i)
+ {
+
+ Body * b = universe.body+i;
+ glColor3f(0.0f, b->mass/1e11*100, 0.0f);
+ //glColor3f(1.0f, 0.0f, 0.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
+ }
+ glFlush();
+
+ AfterDraw(); // Stuff to do after graphics is done drawing
+ // Single-thread - Nothing at all
+ // Multiple threads - signal threads it is safe to change position variables
+
+}
+
+/**
+ * @function Graphics_Keyboard
+ * @purpose This function is to manipulate with the image
+ * @param theKey key pressed
+ * @param mouseX, mouseY position of the mouse in the window
+ */
+void Graphics_Keyboard(unsigned char theKey, int mouseX, int mouseY)
+{
+ if (theKey == 'x' || theKey == 'X')
+ {
+
+ QuitProgram(false);
+
+ return;
+ }
+
+ #ifdef FLYING_CAMERA
+ switch (theKey)
+ {
+ // Translate in direction of camera
+ case 'W':
+ case 'w':
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ eye.p[i] += eye.x[i];
+ break;
+ // Translate backwards from camera direction
+ case 'S':
+ case 's':
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ eye.p[i] -= eye.x[i];
+ break;
+ // Translate left from camera direction
+ case 'A':
+ case 'a':
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ eye.p[i] -= eye.y[i];
+ break;
+ // Translate right from camera direction
+ case 'D':
+ case 'd':
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ eye.p[i] += eye.y[i];
+ break;
+ // Translate up from camera direction
+ case 'Q':
+ case 'q':
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ eye.p[i] += eye.z[i];
+ break;
+ // Translate down from camera direction
+ case 'E':
+ case 'e':
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ eye.p[i] -= eye.z[i];
+ break;
+
+ // Rotate camera direction "down"
+ // (pitch control)
+ case 'I':
+ case 'i':
+ {
+ float theta = M_PI/80.0;
+ Camera old = eye;
+
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ {
+ eye.z[i] = old.z[i] * cos(theta) + old.x[i] * sin(theta);
+ eye.x[i] = old.x[i] * cos(theta) - old.z[i] * sin(theta);
+ }
+ break;
+ }
+ // Rotate camera direction "up"
+ // (pitch control)
+ case 'K':
+ case 'k':
+ {
+ float theta = -M_PI/80.0;
+ Camera old = eye;
+
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ {
+ eye.z[i] = old.z[i] * cos(theta) + old.x[i] * sin(theta);
+ eye.x[i] = old.x[i] * cos(theta) - old.z[i] * sin(theta);
+ }
+ break;
+ }
+ // Rotate camera direction "left"
+ // (yaw control)
+ case 'J':
+ case 'j':
+ {
+ float theta = +M_PI/80.0;
+ Camera old = eye;
+
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ {
+ eye.y[i] = old.y[i] * cos(theta) + old.x[i] * sin(theta);
+ eye.x[i] = old.x[i] * cos(theta) - old.y[i] * sin(theta);
+ }
+ break;
+ }
+ // Rotate camera direction "right"
+ // (yaw control)
+ case 'L':
+ case 'l':
+ {
+ float theta = -M_PI/80.0;
+ Camera old = eye;
+
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ {
+ eye.y[i] = old.y[i] * cos(theta) + old.x[i] * sin(theta);
+ eye.x[i] = old.x[i] * cos(theta) - old.y[i] * sin(theta);
+ }
+ break;
+ }
+ // Rotate camera direction CCW about its axis
+ // (roll control)
+ case 'U':
+ case 'u':
+ {
+ float theta = -M_PI/80.0;
+ Camera old = eye;
+
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ {
+ eye.z[i] = old.z[i] * cos(theta) + old.y[i] * sin(theta);
+ eye.y[i] = old.y[i] * cos(theta) - old.z[i] * sin(theta);
+ }
+ break;
+ }
+ // Rotate camera direction CW about its axis
+ // (roll control)
+ case 'O':
+ case 'o':
+ {
+ float theta = +M_PI/80.0;
+ Camera old = eye;
+
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ {
+ eye.z[i] = old.z[i] * cos(theta) + old.y[i] * sin(theta);
+ eye.y[i] = old.y[i] * cos(theta) - old.z[i] * sin(theta);
+ }
+ break;
+ }
+ }
+
+ /* Code used for debugging the camera
+ system("clear");
+ printf("Camera status:\n");
+ printf("Position: %f %f %f\n", eye.p[0], eye.p[1], eye.p[2]);
+ printf("x: %f %f %f (%f)\n", eye.x[0], eye.x[1], eye.x[2], sqrt(eye.x[0]*eye.x[0] + eye.x[1]*eye.x[1] + eye.x[2]*eye.x[2]));
+ printf("y: %f %f %f (%f)\n", eye.y[0], eye.y[1], eye.y[2], sqrt(eye.y[0]*eye.y[0] + eye.y[1]*eye.y[1] + eye.y[2]*eye.y[2]));
+ printf("z: %f %f %f (%f)\n", eye.z[0], eye.z[1], eye.z[2], sqrt(eye.z[0]*eye.z[0] + eye.z[1]*eye.z[1] + eye.z[2]*eye.z[2]));
+ */
+ #else // The original view code
+
+ // I like how 'I' is used three times.
+
+ if (theKey == 'i' || theKey == 'I') {
+ eyePhi -= M_PI / 20;
+ if (eyePhi == 0)
+ eyePhi = 2 * M_PI;
+ } else if (theKey == 'm' || theKey == 'I') {
+ eyePhi += M_PI / 20;
+ } else if (theKey == 'j' || theKey == 'J') {
+ eyeTheta -= M_PI / 20;
+ } else if (theKey == 'k' || theKey == 'K') {
+ eyeTheta += M_PI / 20;
+ } else if (theKey == ',') {
+ eyeRho += 0.5;
+ } else if (theKey == '.' || theKey == 'I') {
+ eyeRho -= 0.5;
+ } else if (theKey == 'w' || theKey == 'W') {
+ look[1] += 0.5;
+ } else if (theKey == 'z' || theKey == 'Z') {
+ look[1] -= 0.5;
+ } else if (theKey == 'a' || theKey == 'A') {
+ look[0] -= 0.5;
+ } else if (theKey == 's' || theKey == 'S') {
+ look[0] += 0.5;
+ } else if (theKey == '+') {
+ scale *= 1.1;
+ } else if (theKey == '-') {
+ scale *= 0.9;
+ }
+ if (sin(eyePhi) > 0) upY = 1;
+ else upY = 1;
+ #endif //FLYING_CAMERA
+}
+
+/**
+ * @function Graphics_Reshape
+ * @purpose Resize the view window
+ * @param width, height - New size of the window
+ */
+void Graphics_Reshape(int width, int height)
+{
+ double displayRatio = 1.0 * width / height;
+ windowWidth = width;
+ windowHeight = height;
+ glMatrixMode(GL_PROJECTION);
+ glLoadIdentity();
+ gluPerspective(VIEW_ANGLE, displayRatio, WORLD_NEAR, WORLD_FAR);
+ glViewport(0, 0, windowWidth, windowHeight);
+}
+
+
+
#include "../single-thread/nbody.c" // Include original code
#include "graphics.h" // For Graphics_Run function only
-omp_lock_t graphics_lock;
-bool graphics_busy = false;
+
+#ifdef OVER_ENGINEERED
+System * sub_system = NULL;
+#endif //OVER_ENGINEERED
/**
* @function Simulation_Run
*/
void Simulation_Run(int argc, char ** argv)
{
- omp_set_nested(1);
- if (omp_get_nested() != 1)
- {
- fprintf(stderr, "Couldn't set nested parallelism, I kind of need that\n");
- exit(EXIT_FAILURE);
- }
-
-
+
if (options.num_threads == 0)
options.num_threads = omp_get_max_threads();
if (options.draw_graphics)
- omp_set_num_threads(2);
- else
- omp_set_num_threads(1);
-
+ {
+ Graphics_Run(argc, argv);
+ if (options.pedantic_graphics == false)
+ {
+ options.num_threads += 1;
+ omp_set_nested(1);
+ if (omp_get_nested() == 0)
+ {
+ fprintf(stderr, "Couldn't enable nested parallelism\n.");
+ exit(EXIT_FAILURE);
+ }
- omp_init_lock(&graphics_lock);
- #pragma omp parallel sections
- {
-
- #pragma omp section
- {
- omp_set_num_threads(options.num_threads);
- while (true)
+ #pragma omp parallel sections num_threads(2)
{
-
- if (runstate != RUN) break; //Check whether the main thread wants to exit
-
- if (ExitCondition())
+ #pragma omp section
{
- break;
- }
-
- if (options.verbosity != 0 && universe.steps % options.verbosity == 1)
+ Compute();
+ }
+ #pragma omp section
{
- DisplayStatistics();
+ glutMainLoop();
}
-
-
- Compute();
}
- //printf("Left compute loop\n");
+ return;
}
+ }
+ Compute();
- #pragma omp section
- {
- if (options.draw_graphics)
- Graphics_Run(argc, argv);
- //printf("Got to bit after Graphics_Run()\n");
- }
- }
- omp_destroy_lock(&graphics_lock);
}
-/**
- * @function Compute
- * @purpose Compute a single step, multithreaded using OpenMP
- */
-void Compute(void)
+
+void Compute(void)
{
- //printf("There are %d threads working and %d available\n", omp_get_num_threads(), omp_get_max_threads());
- //return;
- unsigned n_per_thread = (universe.N) / options.num_threads;
- unsigned remainder = (universe.N) % options.num_threads;
- #pragma omp parallel for
- for (unsigned i = 0; i < options.num_threads; ++i)
+ omp_set_num_threads(options.num_threads);
+ if (omp_get_max_threads() != options.num_threads)
{
- //printf("Thread %d executes iteration %u\n", omp_get_thread_num(), i);
- System s;
- s.body = universe.body + (i * n_per_thread);
- s.N = (i == options.num_threads - 1) ? n_per_thread : n_per_thread + remainder;
- System_Forces(&s, &universe);
+ fprintf(stderr, "Couldn't use %d threads!\n", options.num_threads);
+ exit(EXIT_FAILURE);
}
- // Wait for graphics to finish drawing
-
- if (options.draw_graphics && options.pedantic_graphics)
- {
- while (graphics_busy);
- }
+ #ifdef OVER_ENGINEERED
+ sub_system = Split_System(&universe, options.num_threads); // I could use omp for loops instead, but I like this way
+ #endif //OVER_ENGINEERED
- #pragma omp parallel for
- for (unsigned i = 0; i < options.num_threads; ++i)
+ #pragma omp parallel
{
- //printf("Thread %d executes iteration %u\n", omp_get_thread_num(), i);
- System s;
- s.body = universe.body + (i * n_per_thread);
- s.N = (i == options.num_threads - 1) ? n_per_thread : n_per_thread + remainder;
- System_Positions(&s);
- }
-
- universe.steps += 1;
+ #ifdef OVER_ENGINEERED
+ int index = omp_get_thread_num();
+ #endif //OVER_ENGINEERED
+ while (!ExitCondition())
+ {
- // Tell graphics to continue
- if (options.draw_graphics && options.pedantic_graphics)
- {
- omp_set_lock(&graphics_lock);
- graphics_busy = true;
- omp_unset_lock(&graphics_lock);
+ #ifdef OVER_ENGINEERED
+ System_Forces(sub_system+index, &universe);
+ #pragma omp barrier // Look! An explicit barrier!
+ System_Positions(sub_system+index);
+ #pragma omp barrier
+ #else
+ #pragma omp for
+ for (unsigned a = 0; a < universe.N; ++a)
+ {
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ (universe.body+a)->F[i] = 0;
+ //printf("Thread %d computes force for body %d\n", omp_get_thread_num(), a);
+ omp_set_num_threads(options.nested_threads);
+ #pragma omp parallel for
+ for (unsigned b = 0; b < universe.N; ++b)
+ {
+ if (b != a)
+ Body_Force(universe.body+a, universe.body+b);
+ }
+ }
- }
+ #pragma omp for
+ for (unsigned a = 0; a < universe.N; ++a)
+ {
+ //printf("Thread %d computes position for body %d\n", omp_get_thread_num(), a);
+ Body_Velocity(universe.body+a);
+ Body_Position(universe.body+a);
+ }
-}
+ #endif //OVER_ENGINEERED
-/**
- * @function BeforeDraw
- * @purpose Called before the graphics thread draws bodies.
- * If pedantic graphics are used, will wait until graphics_busy is set to true before drawing
- */
-void BeforeDraw()
-{
- if (options.verbosity != 0 && universe.steps % options.verbosity == 0)
- DisplayStatistics();
- if (!options.pedantic_graphics)
- return;
- while (graphics_busy == false);
+ #pragma omp master
+ {
+ // The if statement has a cost, but if can be removed easily in an optimised version.
+ // Only the main thread is allowed to mess with glut, because glut is evil (not thread safe).
+ // Although this function may be nested, that will only happen if the if would be false.
+ // ("if the if would be false"; I am awesome at explaining code)
+ if (options.draw_graphics && options.pedantic_graphics)
+ {
+ glutMainLoopEvent();
+ Graphics_Display();
+ }
+ }
+ #pragma omp single
+ {
+ // Only one thread should execute this, and we don't care which one.
+ // So much easier than the pthread version!
+ universe.steps += 1;
+ if (options.verbosity != 0 && universe.steps % options.verbosity == 0)
+ DisplayStatistics();
+ }
+
+ }
+ }
}
-/**
- * @function AfterDraw
- * @purpose Called after the graphics thread draws bodies
- * If pedantic graphics used, will unset graphics_busy when done drawing
- */
-void AfterDraw()
-{
- if (!options.pedantic_graphics)
- return;
- omp_set_lock(&graphics_lock);
- graphics_busy = false;
- omp_unset_lock(&graphics_lock);
-}
--- /dev/null
+/**
+ * @file nbody.c
+ * @purpose 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_Run function only
+
+
+#ifdef OVER_ENGINEERED
+System * sub_system = NULL;
+#endif //OVER_ENGINEERED
+
+/**
+ * @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)
+{
+
+ if (options.num_threads == 0)
+ options.num_threads = omp_get_max_threads();
+
+ if (options.draw_graphics)
+ {
+ Graphics_Run(argc, argv);
+ if (options.pedantic_graphics == false)
+ {
+ options.num_threads += 1;
+ omp_set_nested(1);
+ if (omp_get_nested() == 0)
+ {
+ fprintf(stderr, "Couldn't enable nested parallelism\n.");
+ exit(EXIT_FAILURE);
+ }
+
+
+ #pragma omp parallel sections num_threads(2)
+ {
+ #pragma omp section
+ {
+ Compute();
+ }
+ #pragma omp section
+ {
+ glutMainLoop();
+ }
+ }
+ return;
+ }
+ }
+ #ifdef OVER_ENGINEERED
+ Compute();
+ #else
+ Compute2();
+ #endif //OVER_ENGINEERED
+
+
+
+}
+
+
+void Compute(void)
+{
+ 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);
+ }
+
+ #ifdef OVER_ENGINEERED
+ sub_system = Split_System(&universe, options.num_threads); // I could use omp for loops instead, but I like this way
+ #endif //OVER_ENGINEERED
+
+ #pragma omp parallel
+ {
+ #ifdef OVER_ENGINEERED
+ int index = omp_get_thread_num();
+ #endif //OVER_ENGINEERED
+ while (!ExitCondition())
+ {
+
+ #ifdef OVER_ENGINEERED
+ System_Forces(sub_system+index, &universe);
+ #pragma omp barrier // Look! An explicit barrier!
+ System_Positions(sub_system+index);
+ #pragma omp barrier
+ #else
+ #pragma omp for
+ for (unsigned a = 0; a < universe.N; ++a)
+ {
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ (universe.body+a)->F[i] = 0;
+ //printf("Thread %d computes force for body %d\n", omp_get_thread_num(), a);
+ omp_set_num_threads(options.nested_threads);
+ #pragma omp parallel for
+ for (unsigned b = 0; b < universe.N; ++b)
+ {
+ if (b != a)
+ Body_Force(universe.body+a, universe.body+b);
+ }
+ }
+
+ #pragma omp for
+ for (unsigned a = 0; a < universe.N; ++a)
+ {
+ //printf("Thread %d computes position for body %d\n", omp_get_thread_num(), a);
+ Body_Velocity(universe.body+a);
+ Body_Position(universe.body+a);
+ }
+
+ #endif //OVER_ENGINEERED
+
+ #pragma omp master
+ {
+ // The if statement has a cost, but if can be removed easily in an optimised version.
+ // Only the main thread is allowed to mess with glut, because glut is evil (not thread safe).
+ // Although this function may be nested, that will only happen if the if would be false.
+ // ("if the if would be false"; I am awesome at explaining code)
+ if (options.draw_graphics && options.pedantic_graphics)
+ {
+ glutMainLoopEvent();
+ Graphics_Display();
+ }
+ }
+ #pragma omp single
+ {
+ // Only one thread should execute this, and we don't care which one.
+ // So much easier than the pthread version!
+ universe.steps += 1;
+ if (options.verbosity != 0 && universe.steps % options.verbosity == 0)
+ DisplayStatistics();
+ }
+
+ }
+ }
+}
+
#undef Simulation_Run
void Simulation_Run(int argc, char ** argv);
//#undef QuitProgram
-//void QuitProgram(bool error);
+//#define QuitProgram(x) break
#undef BeforeDraw
-void BeforeDraw();
-#undef AfterDraw
-void AfterDraw();
+#define BeforeDraw() (void)0 // Do nothing (apparently this is how to do nothing with a macro)
+//void BeforeDraw();
-void Compute(); //Compute a single step
+#undef AfterDraw
+#define AfterDraw() (void)0
+//void AfterDraw();
+
+
+
+void Compute(void);
+void Compute2(void);
+
+//#define OVER_ENGINEERED
+// define to (my original approach) manually split the System array (same as pthread version) between threads
+// Manually splitting the System array is not required in openmp, since "#pragma omp for" automatically does it
+// I still like doing it this way, since you can explicitly see what each thread is doing, but it would probably lose me marks.
+// However I can't bring myself to delete it. Hence, compromise.
#endif //_NBODY_OPENMP_H
--- /dev/null
+#ifndef _NBODY_OPENMP_H
+#define _NBODY_OPENMP_H
+
+/**
+ * @file nbody.h
+ * @purpose OpenMP version of N-Body simulator, declarations
+ * @author Sam Moore (20503628) - 2012
+ */
+
+#include "../single-thread/nbody.h" //Include original code
+#include <omp.h>
+
+
+#undef SINGLE_THREADED
+#define OMP_THREADED
+
+
+#define CRAPPY_VERSION
+
+// Replace default macros with thread-safe functions
+#undef Simulation_Run
+void Simulation_Run(int argc, char ** argv);
+//#undef QuitProgram
+//#define QuitProgram(x) break
+
+#undef BeforeDraw
+#define BeforeDraw() (void)0 // Do nothing (apparently this is how to do nothing with a macro)
+//void BeforeDraw();
+
+#undef AfterDraw
+#define AfterDraw() (void)0
+//void AfterDraw();
+
+omp_lock_t graphics_lock;
+
+
+void Compute(void);
+void Compute2(void);
+
+//#define OVER_ENGINEERED
+// define to (my original approach) manually split the System array (same as pthread version) between threads
+// Manually splitting the System array is not required in openmp, since "#pragma omp for" automatically does it
+// I still like doing it this way, since you can explicitly see what each thread is doing, but it would probably lose me marks.
+// However I can't bring myself to delete it. Hence, compromise.
+
+#endif //_NBODY_OPENMP_H
+
+//EOF
//This option must be set, or when glut leaves the main loop. the exit(3) function is called... annoying
glutSetOption(GLUT_ACTION_ON_WINDOW_CLOSE, GLUT_ACTION_CONTINUE_EXECUTION);
- glutMainLoop();
-
+
+ #ifndef OMP_THREADED
+ glutMainLoop();
+ #endif //OMP_THREADED
}
/**
{
if (theKey == 'x' || theKey == 'X')
{
+
QuitProgram(false);
+ glutLeaveMainLoop();
return;
}
* @param a - The Body
* @param b - The System
*/
-inline void Body_Force(Body * a, System * s)
+inline void Body_Forces(Body * a, System * s)
{
- double distance;
- double con;
- double gd;
+
for (unsigned i = 0; i < DIMENSIONS; ++i) //Set Force to zero
a->F[i] = 0;
if (b == a)
continue; //Otherwise we would have infinite force
- //Calculate distance between a and b
- distance = 0.0;
- for (unsigned i = 0; i < DIMENSIONS; ++i)
- distance += square(b->x[i] - a->x[i]);
- distance = sqrt(distance);
- con = G * a->mass * b->mass / square(distance);
- gd = con / distance;
- for (unsigned i = 0; i < DIMENSIONS; ++i) //Add force from b to a's net force
- a->F[i] += gd * (b->x[i] - a->x[i]);
+ Body_Force(a, b);
}
}
+/**
+ * @function Body_Force
+ * @purpose Compute force on body a due to body b
+ * @param a - Body a
+ * @param b - Body b
+ */
+inline void Body_Force(Body * a, Body * b)
+{
+ double distance;
+ double con;
+ double gd;
+ //Calculate distance between a and b
+ distance = 0.0;
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ distance += square(b->x[i] - a->x[i]);
+ distance = sqrt(distance);
+ con = G * a->mass * b->mass / square(distance);
+ gd = con / distance;
+ for (unsigned i = 0; i < DIMENSIONS; ++i) //Add force from b to a's net force
+ a->F[i] += gd * (b->x[i] - a->x[i]);
+}
+
/**
* @function Body_Velocity
* @purpose Compute velocity of a body
//printf("Compute forces for %p - %d bodies...\n", (void*)s1, s1->N);
for (unsigned i = 0; i < s1->N; ++i)
{
- Body_Force(s1->body+i, s2);
+ Body_Forces(s1->body+i, s2);
Body_Velocity(s1->body+i);
}
//printf("Compute positions for %p - Done!\n", (void*)s1);
//fprintf(stderr,"runstate %d\n timeout %d\n steps %d\n", (int)(runstate != RUN), (int)(options.timeout > 0.00 && ((unsigned)(time(NULL) - options.start_time.tv_sec) >= options.timeout)), (int)(options.num_steps > 0 && universe.steps > options.num_steps));
return result;
}
+
+/**
+ * @function Split_System
+ * @purpose Helper to divide one system into an array of systems
+ * Each sub system will have N = (s->N / n) bodies in it
+ * @param s - The original system (typically &universe)
+ * @param n - The number of sub systems in the array
+ *
+ * WARNING: It is the caller's responsibility to free() the returned array
+ */
+System * Split_System(System * s, unsigned n)
+{
+ System * result = (System*)(calloc(n, sizeof(System)));
+ if (result == NULL)
+ {
+ perror("Couldn't create array of sub systems");
+ exit(EXIT_FAILURE);
+ }
+
+ unsigned n_per_system = (s->N) / n;
+ unsigned remainder = (s->N) % n;
+
+ for (unsigned i = 0; i < n; ++i)
+ {
+ result[i].N = n_per_system;
+ if (i == n-1)
+ result[i].N += remainder;
+ result[i].body = (s->body) + (n_per_system * i);
+ result[i].steps = 0;
+ }
+ return result;
+}
--- /dev/null
+/**
+ * @file nbody.c
+ * @author Sam Moore (20503628) 2012
+ * @purpose N-Body simulator - Definition of simulation functions; single threaded version
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <time.h>
+#include <string.h>
+#include <stdbool.h>
+#include <stdlib.h>
+#include <unistd.h>
+
+#include "nbody.h"
+
+RUNSTATE runstate = RUN;
+
+/**
+ * @function Body_Print
+ * @purpose Print the body to a file (or stdout)
+ * @param a - Body to print'
+ * @param out - FILE* to print to
+ */
+inline void Body_Print(Body * a, FILE * out)
+{
+ //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",
+ //(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]);
+ fprintf(out, "%f %f %f %f %f %f %f %f %f %f\n",
+ a->mass, a->x[0], a->x[1], a->x[2], a->v[0], a->v[1], a->v[2], a->F[0], a->F[1], a->F[2]);
+}
+
+/**
+ * @function Body_Force
+ * @purpose Calculates the force on a single Body due to all bodies in a System
+ * @param a - The Body
+ * @param b - The System
+ */
+inline void Body_Forces(Body * a, System * s)
+{
+
+
+ for (unsigned i = 0; i < DIMENSIONS; ++i) //Set Force to zero
+ a->F[i] = 0;
+
+ for (unsigned index = 0; index < s->N; ++index) // Step through all bodies in the System
+ {
+ Body * b = s->body+index; //Current body
+ if (b == a)
+ continue; //Otherwise we would have infinite force
+
+ Body_Force(a, b);
+ }
+}
+
+/**
+ * @function Body_Force
+ * @purpose Compute force on body a due to body b
+ * @param a - Body a
+ * @param b - Body b
+ */
+inline void Body_Force(Body * a, Body * b)
+{
+ //Calculate distance between a and b
+ distance = 0.0;
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ distance += square(b->x[i] - a->x[i]);
+ distance = sqrt(distance);
+ con = G * a->mass * b->mass / square(distance);
+ gd = con / distance;
+ for (unsigned i = 0; i < DIMENSIONS; ++i) //Add force from b to a's net force
+ a->F[i] += gd * (b->x[i] - a->x[i]);
+}
+
+/**
+ * @function Body_Velocity
+ * @purpose Compute velocity of a body
+ * @param a - The Body
+ */
+inline void Body_Velocity(Body * a)
+{
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ a->v[i] += a->F[i] / a->mass * DELTA_T;
+}
+
+/**
+ * @function Body_Position
+ * @purpose Compute position of a body
+ * @param a - The Body
+ */
+inline void Body_Position(Body * a)
+{
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ a->x[i] += a->v[i] * DELTA_T;
+}
+
+/**
+ * @function System_Compute
+ * @purpose Compute forces and then positions for bodies in System
+ * NOTE: Only used in the single threaded version of the program
+ */
+inline void System_Compute(System * s)
+{
+ System_Forces(s, s);
+ System_Positions(s);
+}
+
+/**
+ * @function System_Forces
+ * @purpose Compute forces on each object in System s1 due to all bodies in System s2
+ * @param s1, s2 - The two systems
+ * NOTE: In single threaded version, called with s1 == s2 == &universe
+ * In multi threaded version, called with s2 == &universe, but s1 is constructed for each thread
+ * In nested multi-threaded version, s2 is constructed for the nested threads as well
+ */
+inline void System_Forces(System * s1, System * s2)
+{
+ //printf("Compute forces for %p - %d bodies...\n", (void*)s1, s1->N);
+ for (unsigned i = 0; i < s1->N; ++i)
+ {
+ Body_Forces(s1->body+i, s2);
+ Body_Velocity(s1->body+i);
+ }
+ //printf("Compute positions for %p - Done!\n", (void*)s1);
+}
+
+/**
+ * @function System_Positions
+ * @purpose Compute positions of all objects in System
+ * @param s - The system
+ */
+inline void System_Positions(System * s)
+{
+ //printf("Compute positions for %p - %d bodies...\n", (void*)s, s->N);
+ for (unsigned i = 0; i < s->N; ++i)
+ Body_Position(s->body+i);
+ //printf("Compute positions for %p - Done!\n", (void*)s);
+
+}
+
+
+
+
+/**
+ * @function System_Init
+ * @purpose Initialise a System from an input file
+ * @param s - The System
+ * @param fileName - The input file
+ */
+inline void System_Init(System * s, const char *fileName)
+{
+ char line[LINE_SIZE];
+ char * token;
+ FILE * file;
+
+ file = fopen(fileName, "rt");
+ s->N = atoi(fgets(line, LINE_SIZE, file));
+ s->body = (Body*) calloc((size_t)s->N, sizeof(Body));
+ s->steps = 0;
+
+ //printf("----------------------Initial field-------------------------------\n");
+
+ for (unsigned i = 0; i < s->N; ++i)
+ {
+ if (fgets(line, LINE_SIZE, file) != NULL)
+ {
+ Body * a = s->body+i;
+ token = strtok(line, ",; \t");
+ a->mass = atof(token);
+
+ for (unsigned j = 0; j < DIMENSIONS; ++j)
+ {
+ token = strtok(NULL, ",; \t");
+ a->x[j] = atof(token);
+ }
+ for (unsigned j = 0; j < DIMENSIONS; ++j)
+ {
+ token = strtok(NULL, ",; \t");
+ a->v[j] = atof(token);
+ }
+
+ // Ignore any additional tokens
+ while (token != NULL)
+ token = strtok(NULL, ",; \t");
+ //Body_Print(a);
+ }
+ }
+
+ //printf("--------------------------------------------------------------\n");
+
+ fclose(file);
+}
+
+/**
+ * @function Universe_Cleanup
+ * @purpose Cleans up the universe by freeing all memory
+ * Also prints the bodies in the universe to a file specified in "options" if required.
+ */
+inline void Universe_Cleanup()
+{
+ //fprintf(stderr, "Called Universe_Cleanup()\n");
+ if (options.output != NULL) // An output file was specified...
+ {
+ FILE * save = NULL;
+ if (strcmp("stdout", options.output) == 0)
+ save = stdout;
+ else if (strcmp("stderr", options.output) == 0)
+ save = stderr;
+ else
+ save = fopen(options.output, "w");
+ if (save == NULL)
+ perror("Couldn't save universe to file.");
+ else
+ {
+ // Print the output in the same format as the initial field file
+ fprintf(save, "%u\n", universe.N);
+
+ for (unsigned i = 0; i < universe.N; ++i) // Print all bodies to the file
+ {
+ Body_Print(universe.body+i, save);
+
+
+ }
+ fclose(save);
+ }
+
+ }
+ free(universe.body); // Cleanup memory used for bodies
+
+}
+
+/**
+ * @function DisplayStatistics
+ * @purpose Display useful information about the computations done,
+ * - Total number of steps computed
+ * - Time in seconds elapsed
+ * - Clock cycles executed
+ * - Equivelant time for a single thread to execute same number of clock cycles
+ */
+inline void DisplayStatistics()
+{
+ clock_t end = clock();
+ struct timeval end_time;
+ if (gettimeofday(&end_time, NULL) != 0)
+ {
+ perror("Couldn't get time of day.\n");
+ return;
+ }
+ float runtime = (float)(end_time.tv_sec-options.start_time.tv_sec) + (float)(end_time.tv_usec-options.start_time.tv_usec) / 1.0e6;
+ float clock_time = (float)(end - options.start_clock) / (float)(CLOCKS_PER_SEC);
+ printf("%u\t%f\t%u\t%f\n", universe.steps, runtime, (unsigned)(end - options.start_clock), clock_time);
+}
+
+
+/**
+ * @function ExitCondition
+ * @purpose Helper to check whether the program is supposed to exit
+ * Does not check for user triggered quit
+ * Checks for either a timeout, or computation of the required number of steps
+ *
+ * The reason timeouts are integers is because this function is called a lot in lots of threads
+ * So using floats (and working out microseconds every time) is expensive
+ *
+ */
+inline bool ExitCondition(void)
+{
+ bool result = (runstate != RUN ||
+ (options.num_steps >= 0 && universe.steps >= options.num_steps) ||
+ (options.timeout >= 0 && ((int)(time(NULL) - options.start_time.tv_sec) >= options.timeout)));
+ //fprintf(stderr,"runstate %d\n timeout %d\n steps %d\n", (int)(runstate != RUN), (int)(options.timeout > 0.00 && ((unsigned)(time(NULL) - options.start_time.tv_sec) >= options.timeout)), (int)(options.num_steps > 0 && universe.steps > options.num_steps));
+ return result;
+}
+
+/**
+ * @function Split_System
+ * @purpose Helper to divide one system into an array of systems
+ * Each sub system will have N = (s->N / n) bodies in it
+ * @param s - The original system (typically &universe)
+ * @param n - The number of sub systems in the array
+ *
+ * WARNING: It is the caller's responsibility to free() the returned array
+ */
+System * Split_System(System * s, unsigned n)
+{
+ System * result = (System*)(calloc(n, sizeof(System)));
+ if (result == NULL)
+ {
+ perror("Couldn't create array of sub systems");
+ exit(EXIT_FAILURE);
+ }
+
+ unsigned n_per_system = (s->N) / n;
+ unsigned remainder = (s->N) % n;
+
+ for (unsigned i = 0; i < n; ++i)
+ {
+ result[i].N = n_per_system;
+ if (i == n-1)
+ result[i].N += remainder;
+ result[i].body = (s->body) + (n_per_system * i);
+ result[i].steps = 0;
+ }
+ return result;
+}
} Options;
void Body_Print(Body * a, FILE * out); //Print body a
-void Body_Force(Body * a, System * s); //Compute force on body a due to system of bodies s
+void Body_Forces(Body * a, System * s); //Compute force on body a due to system of bodies s
+void Body_Force(Body * a, Body * b); // Compute force on body a due to body b
void Body_Velocity(Body * a); //Compute velocity of body a
void Body_Position(Body * a); //Compute position of body a
void System_Compute(System * s);
void System_Forces(System * s1, System * s2); //Compute forces for bodies in s1 due to bodies in s2 (also updates velocities)
void System_Positions(System * s); //Update positions for bodies in s1
-
+System * Split_System(System * s, unsigned n); // Splits one system into a number of other systems, returns an array of size n
void Universe_Cleanup(); //Cleanup universe and write bodies to file
--- /dev/null
+#ifndef _NBODY_H
+#define _NBODY_H
+
+/**
+ * @file nbody.h
+ * @author Sam Moore (205030628)
+ * @purpose N-Body simulator: declarations of simulation related parameters
+ */
+
+#include <stdlib.h>
+#include <stdbool.h>
+#include <stdio.h>
+#include <time.h>
+#include <sys/time.h> //POSIX time
+
+
+
+#define SINGLE_THREADED
+
+// --- The below macros will be undefined by the multithreaded versions, and replaced with functions --- //
+
+//Sets up the simulation; in multithreaded versions, will spawn threads
+#define Simulation_Run(argc, argv) \
+if (options.draw_graphics) \
+ Graphics_Run(argc, argv); \
+else \
+{ \
+ while (ExitCondition() == false) \
+ { BeforeDraw(); AfterDraw(); } \
+}
+
+#define QuitProgram(error) (runstate = (error == true) ? QUIT : QUIT_ERROR) //Prepares to exit program, is thread safe in multithreaded versions
+
+//Macro to be overwritten in multithreaded versions, called before the graphics is allowed to draw anything
+#define BeforeDraw() \
+System_Compute(&universe); \
+universe.steps += 1; \
+if (options.verbosity != 0 && universe.steps % options.verbosity == 0) \
+ DisplayStatistics(); \
+
+
+
+
+//Macro to be overwritten in multithreaded versions, called after the graphics has finished drawing.
+#define AfterDraw()
+
+// --- Constants and other Macros --- //
+
+#define M_PI 3.14159265358979323846264338327950288 /* pi */
+#define G 6.67428E-11
+#define DELTA_T 0.05
+#define DIMENSIONS 3
+#define LINE_SIZE 1000
+#define square(x) ((x)*(x))
+
+
+/**
+ * Structure to represent a single Body
+ * @param mass - Mass of the body
+ * @param x - Position vector (array)
+ * @param v - Velocity vector (array)
+ * @param F - Net force vector (array)
+ */
+typedef struct
+{
+
+ double mass;
+ double x[DIMENSIONS];
+ double v[DIMENSIONS];
+ double F[DIMENSIONS];
+
+} Body;
+
+/**
+ * Structure to store an array of bodies, along with the size of the array.
+ * The universe is represented by a single System.
+ * In the multithreaded program, the universe is subdivided into one system for each working thread to use.
+ * @param N - Size of the array
+ * @param body - The array of bodies
+ */
+typedef struct
+{
+ unsigned N; // Number of bodies in the System
+ Body * body; // Array of bodies
+
+ unsigned steps; //Number of steps simulated
+
+} System;
+
+/**
+ * Structure to represent options passed to the program.
+ */
+typedef struct
+{
+ const char * input; // initial body field
+ const char * output; // file to write final positions / velocities of bodies to
+ const char * program; // program name
+ int num_threads; // number of worker threads to spawn (must be greater than 1 for any to be spawned)
+ int nested_threads; // number of threads to nest computations with (must be greater than 1 for any to be spawned)
+ int num_steps; // number of steps to run before stopping (run indefinately if less than zero)
+ int timeout; // number of seconds to run before stopping (run indefinately if less than zero)
+ bool draw_graphics; // whether or not to actually draw graphics
+ bool pedantic_graphics; // whether the graphics thread will synchronise with the computation thread (true) or just draw as fast as possible (false)
+ bool print_positions; // print positions of bodies to stdout on every step
+ int verbosity; // print statistics every number of steps indicated by this variable
+ clock_t start_clock; // clock cycles done when simulation starts
+ struct timeval start_time; // time at which simulation starts
+} Options;
+
+void Body_Print(Body * a, FILE * out); //Print body a
+void Body_Force(Body * a, System * s); //Compute force on body a due to system of bodies s
+void Body_Force(Body * a, Body * b); // Compute force on body a due to body b
+void Body_Velocity(Body * a); //Compute velocity of body a
+void Body_Position(Body * a); //Compute position of body a
+
+void System_Init(System * s, const char * fileName); //Initialise System (array of bodies) from a text file
+void System_Compute(System * s);
+void System_Forces(System * s1, System * s2); //Compute forces for bodies in s1 due to bodies in s2 (also updates velocities)
+void System_Positions(System * s); //Update positions for bodies in s1
+System * Split_System(System * s, unsigned n); // Splits one system into a number of other systems, returns an array of size n
+
+void Universe_Cleanup(); //Cleanup universe and write bodies to file
+
+void DisplayStatistics(); // Print information about steps computed, total (real) runtime, and CPU cycles
+
+
+bool ExitCondition(void); //Checks whether the program is supposed to exit automatically yet (ie: other than the user pressing "quit")
+
+
+typedef enum {RUN, QUIT, QUIT_ERROR} RUNSTATE;
+extern RUNSTATE runstate; // Set runstate to QUIT or QUIT_ERROR and the simulation will stop on the next step.
+// This is fairly redundant in the single threaded version, but useful for making sure *all* threads know to exit in the multi-threaded version
+
+
+extern System universe; // The main array of bodies; global variable.
+extern Options options; // Parameters passed to program
+
+
+
+#endif //_NBODY_H
git.ucc.asn.au / matches / honours.git / commitdiff
UCC git Repository :: git.ucc.asn.au