CPSC 521 assignment 2 source
Go back. Source for
nbody.c from assignment 2:
/* This program is written in C99. */
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h> // for memcpy
#include <math.h> // for hypot*
#define USAGE \
"Usage: %s rounds granularity inputfile\n\n" \
"This is an n-body simulation program, parallelized with MPI.\n"
#define SCALAR double
#define SCALAR_MPI_TYPE MPI_DOUBLE
#if SCALAR == float
#define SQRT sqrtf
#define HYPOT hypotf
#else
#define SQRT sqrt
#define HYPOT hypot
#endif
//#define DEBUG_OUTPUT
//#define VISUALIZE_OUTPUT
#define DO_TIMING
//#define DO_TIMING_FINE
#ifdef DEBUG_OUTPUT
#define DEBUG(s, ...) printf(s, __VA_ARGS__)
#else
#define DEBUG(s, ...)
#endif
enum mpi_tag_t {
TAG_INIT,
TAG_SIMULATE,
TAG_TERMINATE
};
typedef struct {
SCALAR x, y;
} vector_t;
typedef struct {
vector_t pos, velocity;
SCALAR mass;
} body_t;
body_t *parse_input(const char *filename, int lines);
void setup_simulation(body_t *data, int rounds, int granularity,
int self, int bodies);
void perform_simulation(body_t *body_array, int rounds, int granularity,
int self, int bodies);
void simulation_step(body_t *body_array, int granularity, int self, int bodies
#ifdef VISUALIZE_OUTPUT
, int visualize_round
#endif
);
void send_body_array(body_t *body, int count, int dest, int tag);
void recv_body_array(body_t *body, int count, int source, int tag);
int main(int argc, char *argv[]) {
MPI_Init(&argc, &argv);
#ifdef DO_TIMING
double start_time = MPI_Wtime();
#endif
int rank, size;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
if(argc != 4) {
if(rank == 0) printf(USAGE, argv[0]);
return 1;
}
else {
int rounds = strtol(argv[1], NULL, 0);
int granularity = strtol(argv[2], NULL, 0);
body_t *data = NULL;
if(rank == 0) {
data = parse_input(argv[3], size*granularity);
}
setup_simulation(data, rounds, granularity, rank, size);
free(data); // do nothing if NULL
}
#ifdef DO_TIMING
if(rank == 0) {
fprintf(stderr, "Total wall time %f\n", MPI_Wtime() - start_time);
}
#endif
MPI_Finalize();
return 0;
}
body_t *parse_input(const char *filename, int lines) {
FILE *fp = fopen(filename, "rt");
if(!fp) {
printf("Can't open input file \"%s\"\n", filename);
exit(1);
}
body_t *data = malloc(lines * sizeof(*data));
if(!data) {
printf("Can't allocate memory for %d bodies\n", lines);
fclose(fp);
exit(1);
}
char buffer[BUFSIZ];
int i;
for(i = 0; i < lines && fgets(buffer, sizeof buffer, fp); i ++) {
double x, y, mass;
if(sscanf(buffer, "%lf%lf%lf", &x, &y, &mass) != 3) {
printf("Parse error on input line %d\n", i + 1);
exit(1);
}
data[i].pos.x = x;
data[i].pos.y = y;
data[i].velocity.x = 0.0;
data[i].velocity.y = 0.0;
data[i].mass = mass;
}
fclose(fp);
if(i != lines) {
printf("Not enough bodies in input: found %d, need %d\n", i, lines);
exit(1);
}
return data;
}
void setup_simulation(body_t *data, int rounds, int granularity,
int self, int bodies) {
body_t *body_array = malloc(granularity * sizeof(*body_array));
if(self == 0) {
for(int i = 1; i < bodies; i ++) {
send_body_array(data + i*granularity, granularity,
self + 1, TAG_INIT);
}
memcpy(body_array, data + 0, granularity * sizeof(*body_array));
}
else {
recv_body_array(body_array, granularity, self - 1, TAG_INIT);
body_t *temp_body = malloc(granularity * sizeof(*temp_body));
for(int i = self + 1; i < bodies; i ++) {
recv_body_array(temp_body, granularity, self - 1, TAG_INIT);
send_body_array(temp_body, granularity, self + 1, TAG_INIT);
}
free(temp_body);
}
perform_simulation(body_array, rounds, granularity, self, bodies);
if(self == 0) {
for(int j = 0; j < granularity; j ++) {
printf("%.10f %.10f %f\n",
body_array[j].pos.x, body_array[j].pos.y, body_array[j].mass);
}
body_t *temp_body = malloc(granularity * sizeof(*temp_body));
for(int i = 1; i < bodies; i ++) {
recv_body_array(temp_body, granularity, 1, TAG_TERMINATE);
for(int j = 0; j < granularity; j ++) {
printf("%.10f %.10f %f\n",
temp_body[j].pos.x, temp_body[j].pos.y, temp_body[j].mass);
}
}
free(temp_body);
}
else {
send_body_array(body_array, granularity, self - 1, TAG_TERMINATE);
body_t *temp_body = malloc(granularity * sizeof(*temp_body));
for(int i = self + 1; i < bodies; i ++) {
recv_body_array(temp_body, granularity, self + 1, TAG_TERMINATE);
send_body_array(temp_body, granularity, self - 1, TAG_TERMINATE);
}
free(temp_body);
}
free(body_array);
}
void perform_simulation(body_t *body_array, int rounds, int granularity,
int self, int bodies) {
#ifdef VISUALIZE_OUTPUT
if(self == 0) printf("%d\n", bodies*granularity);
#endif
int i;
for(i = 0; i < rounds; i ++) {
if(self == 0) {
DEBUG("--- begin round %d/%d\n", i, rounds);
}
#ifdef VISUALIZE_OUTPUT
if(rounds < 1000 || (i > 0 && i % (rounds / 1000) == 0)) {
simulation_step(body_array, granularity, self, bodies, self == 0);
}
else {
simulation_step(body_array, granularity, self, bodies, 0);
}
#else
simulation_step(body_array, granularity, self, bodies);
#endif
}
}
void simulation_step(body_t *body_array, int granularity, int self, int bodies
#ifdef VISUALIZE_OUTPUT
, int visualize_round
#endif
) {
SCALAR *force_x = malloc(granularity * sizeof(*force_x));
SCALAR *force_y = malloc(granularity * sizeof(*force_y));
const SCALAR TIME_DELTA = 1.0;
const SCALAR CONSTANT_G = 6.6738480e-11; // from empirical physics
for(int a = 0; a < granularity; a ++) {
force_x[a] = force_y[a] = 0.0;
}
#ifdef VISUALIZE_OUTPUT
if(visualize_round) for(int v = 0; v < granularity; v ++) {
printf("%f %f %f\n", body_array[v].pos.x, body_array[v].pos.y, body_array[v].mass);
}
#endif
int prev = (self + bodies - 1) % bodies;
int next = (self + 1) % bodies;
int p;
for(p = 0; p < bodies; p ++) {
if(p == 0) {
if(bodies > 1) send_body_array(body_array, granularity, next, TAG_SIMULATE);
continue;
}
body_t *other = malloc(granularity * sizeof(*other));
recv_body_array(other, granularity, prev, TAG_SIMULATE);
#ifdef VISUALIZE_OUTPUT
if(visualize_round) for(int v = 0; v < granularity; v ++) {
printf("%f %f %f\n", other[v].pos.x, other[v].pos.y, other[v].mass);
}
#endif
if(p + 1 < bodies) send_body_array(other, granularity, next, TAG_SIMULATE);
for(int a = 0; a < granularity; a ++) {
for(int b = 0; b < granularity; b ++) {
vector_t pos_vector;
pos_vector.x = other[b].pos.x - body_array[a].pos.x;
pos_vector.y = other[b].pos.y - body_array[a].pos.y;
SCALAR length = HYPOT(pos_vector.x, pos_vector.y);
if(length > 0.0) {
SCALAR force = (CONSTANT_G * body_array[a].mass * other[b].mass)
/ (length * length);
force_x[a] += force * (pos_vector.x / length);
force_y[a] += force * (pos_vector.y / length);
}
}
}
free(other);
}
for(int a = 0; a < granularity; a ++) {
for(int b = a + 1; b < granularity; b ++) {
vector_t pos_vector;
pos_vector.x = body_array[b].pos.x - body_array[a].pos.x;
pos_vector.y = body_array[b].pos.y - body_array[a].pos.y;
SCALAR length = HYPOT(pos_vector.x, pos_vector.y);
if(length > 0.0) {
SCALAR force
= (CONSTANT_G * body_array[a].mass * body_array[b].mass)
/ (length * length);
force_x[a] += force * (pos_vector.x / length);
force_y[a] += force * (pos_vector.y / length);
force_x[b] += force * (pos_vector.x / length);
force_y[b] += force * (pos_vector.y / length);
}
}
}
for(int a = 0; a < granularity; a ++) {
body_array[a].velocity.x += force_x[a] * TIME_DELTA / body_array[a].mass;
body_array[a].velocity.y += force_y[a] * TIME_DELTA / body_array[a].mass;
body_array[a].pos.x += body_array[a].velocity.x * TIME_DELTA;
body_array[a].pos.y += body_array[a].velocity.y * TIME_DELTA;
DEBUG("node %d body %d is (%f,%f) (%f,%f) %f\n", self, a,
body_array[a].pos.x, body_array[a].pos.y,
body_array[a].velocity.x, body_array[a].velocity.y, body_array[a].mass);
}
free(force_x);
free(force_y);
}
void send_body_array(body_t *body, int count, int dest, int tag) {
SCALAR *array = malloc(count * 3 * sizeof(*array));
for(int i = 0; i < count; i ++) {
array[3*i + 0] = body[i].pos.x;
array[3*i + 1] = body[i].pos.y;
array[3*i + 2] = body[i].mass;
}
MPI_Send(array, count * 3, SCALAR_MPI_TYPE,
dest, tag, MPI_COMM_WORLD);
free(array);
}
void recv_body_array(body_t *body, int count, int source, int tag) {
SCALAR *array = malloc(count * 3 * sizeof(*array));
MPI_Status status;
MPI_Recv(array, count * 3, SCALAR_MPI_TYPE,
source, tag, MPI_COMM_WORLD, &status);
for(int i = 0; i < count; i ++) {
body[i].pos.x = array[3*i + 0];
body[i].pos.y = array[3*i + 1];
body[i].mass = array[3*i + 2];
body[i].velocity.x = 0.0;
body[i].velocity.y = 0.0;
}
free(array);
}