3 * @author Sam Moore (20503628) 2012
4 * @purpose N-Body simulator - Definition of simulation functions; single threaded version
18 RUNSTATE runstate = RUN;
21 * @function Body_Print
22 * @purpose Print the body to a file (or stdout)
23 * @param a - Body to print'
24 * @param out - FILE* to print to
26 inline void Body_Print(Body * a, FILE * out)
28 //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",
29 //(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]);
30 fprintf(out, "%f %f %f %f %f %f %f %f %f %f\n",
31 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]);
35 * @function Body_Force
36 * @purpose Calculates the force on a single Body due to all bodies in a System
38 * @param b - The System
40 inline void Body_Forces(Body * a, System * s)
44 for (unsigned i = 0; i < DIMENSIONS; ++i) //Set Force to zero
47 for (unsigned index = 0; index < s->N; ++index) // Step through all bodies in the System
49 Body * b = s->body+index; //Current body
51 continue; //Otherwise we would have infinite force
58 * @function Body_Force
59 * @purpose Compute force on body a due to body b
63 inline void Body_Force(Body * a, Body * b)
65 //Calculate distance between a and b
67 for (unsigned i = 0; i < DIMENSIONS; ++i)
68 distance += square(b->x[i] - a->x[i]);
69 distance = sqrt(distance);
70 con = G * a->mass * b->mass / square(distance);
72 for (unsigned i = 0; i < DIMENSIONS; ++i) //Add force from b to a's net force
73 a->F[i] += gd * (b->x[i] - a->x[i]);
77 * @function Body_Velocity
78 * @purpose Compute velocity of a body
81 inline void Body_Velocity(Body * a)
83 for (unsigned i = 0; i < DIMENSIONS; ++i)
84 a->v[i] += a->F[i] / a->mass * DELTA_T;
88 * @function Body_Position
89 * @purpose Compute position of a body
92 inline void Body_Position(Body * a)
94 for (unsigned i = 0; i < DIMENSIONS; ++i)
95 a->x[i] += a->v[i] * DELTA_T;
99 * @function System_Compute
100 * @purpose Compute forces and then positions for bodies in System
101 * NOTE: Only used in the single threaded version of the program
103 inline void System_Compute(System * s)
110 * @function System_Forces
111 * @purpose Compute forces on each object in System s1 due to all bodies in System s2
112 * @param s1, s2 - The two systems
113 * NOTE: In single threaded version, called with s1 == s2 == &universe
114 * In multi threaded version, called with s2 == &universe, but s1 is constructed for each thread
115 * In nested multi-threaded version, s2 is constructed for the nested threads as well
117 inline void System_Forces(System * s1, System * s2)
119 //printf("Compute forces for %p - %d bodies...\n", (void*)s1, s1->N);
120 for (unsigned i = 0; i < s1->N; ++i)
122 Body_Forces(s1->body+i, s2);
123 Body_Velocity(s1->body+i);
125 //printf("Compute positions for %p - Done!\n", (void*)s1);
129 * @function System_Positions
130 * @purpose Compute positions of all objects in System
131 * @param s - The system
133 inline void System_Positions(System * s)
135 //printf("Compute positions for %p - %d bodies...\n", (void*)s, s->N);
136 for (unsigned i = 0; i < s->N; ++i)
137 Body_Position(s->body+i);
138 //printf("Compute positions for %p - Done!\n", (void*)s);
146 * @function System_Init
147 * @purpose Initialise a System from an input file
148 * @param s - The System
149 * @param fileName - The input file
151 inline void System_Init(System * s, const char *fileName)
153 char line[LINE_SIZE];
157 file = fopen(fileName, "rt");
158 s->N = atoi(fgets(line, LINE_SIZE, file));
159 s->body = (Body*) calloc((size_t)s->N, sizeof(Body));
162 //printf("----------------------Initial field-------------------------------\n");
164 for (unsigned i = 0; i < s->N; ++i)
166 if (fgets(line, LINE_SIZE, file) != NULL)
168 Body * a = s->body+i;
169 token = strtok(line, ",; \t");
170 a->mass = atof(token);
172 for (unsigned j = 0; j < DIMENSIONS; ++j)
174 token = strtok(NULL, ",; \t");
175 a->x[j] = atof(token);
177 for (unsigned j = 0; j < DIMENSIONS; ++j)
179 token = strtok(NULL, ",; \t");
180 a->v[j] = atof(token);
183 // Ignore any additional tokens
184 while (token != NULL)
185 token = strtok(NULL, ",; \t");
190 //printf("--------------------------------------------------------------\n");
196 * @function Universe_Cleanup
197 * @purpose Cleans up the universe by freeing all memory
198 * Also prints the bodies in the universe to a file specified in "options" if required.
200 inline void Universe_Cleanup()
202 //fprintf(stderr, "Called Universe_Cleanup()\n");
203 if (options.output != NULL) // An output file was specified...
206 if (strcmp("stdout", options.output) == 0)
208 else if (strcmp("stderr", options.output) == 0)
211 save = fopen(options.output, "w");
213 perror("Couldn't save universe to file.");
216 // Print the output in the same format as the initial field file
217 fprintf(save, "%u\n", universe.N);
219 for (unsigned i = 0; i < universe.N; ++i) // Print all bodies to the file
221 Body_Print(universe.body+i, save);
229 free(universe.body); // Cleanup memory used for bodies
234 * @function DisplayStatistics
235 * @purpose Display useful information about the computations done,
236 * - Total number of steps computed
237 * - Time in seconds elapsed
238 * - Clock cycles executed
239 * - Equivelant time for a single thread to execute same number of clock cycles
241 inline void DisplayStatistics()
243 clock_t end = clock();
244 struct timeval end_time;
245 if (gettimeofday(&end_time, NULL) != 0)
247 perror("Couldn't get time of day.\n");
250 float runtime = (float)(end_time.tv_sec-options.start_time.tv_sec) + (float)(end_time.tv_usec-options.start_time.tv_usec) / 1.0e6;
251 float clock_time = (float)(end - options.start_clock) / (float)(CLOCKS_PER_SEC);
252 printf("%u\t%f\t%u\t%f\n", universe.steps, runtime, (unsigned)(end - options.start_clock), clock_time);
257 * @function ExitCondition
258 * @purpose Helper to check whether the program is supposed to exit
259 * Does not check for user triggered quit
260 * Checks for either a timeout, or computation of the required number of steps
262 * The reason timeouts are integers is because this function is called a lot in lots of threads
263 * So using floats (and working out microseconds every time) is expensive
266 inline bool ExitCondition(void)
268 bool result = (runstate != RUN ||
269 (options.num_steps >= 0 && universe.steps >= options.num_steps) ||
270 (options.timeout >= 0 && ((int)(time(NULL) - options.start_time.tv_sec) >= options.timeout)));
271 //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));
276 * @function Split_System
277 * @purpose Helper to divide one system into an array of systems
278 * Each sub system will have N = (s->N / n) bodies in it
279 * @param s - The original system (typically &universe)
280 * @param n - The number of sub systems in the array
282 * WARNING: It is the caller's responsibility to free() the returned array
284 System * Split_System(System * s, unsigned n)
286 System * result = (System*)(calloc(n, sizeof(System)));
289 perror("Couldn't create array of sub systems");
293 unsigned n_per_system = (s->N) / n;
294 unsigned remainder = (s->N) % n;
296 for (unsigned i = 0; i < n; ++i)
298 result[i].N = n_per_system;
300 result[i].N += remainder;
301 result[i].body = (s->body) + (n_per_system * i);