CPSC 521 assignment 1 source

Go back. Source for nbody.c from assignment 1:
/* This program is written in C99. */

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>  // for hypot*

#define USAGE \
    "Usage: %s rounds 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 SERIALIZE_VELOCITIES
//#define VISUALIZE_OUTPUT

#define RING_COMM

#ifdef RING_COMM
    #define RING_COMM_INIT
    #define RING_COMM_SIMULATE
    #define RING_COMM_TERMINATE
#endif

#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 self, int bodies);
void perform_simulation(body_t *body, int rounds, int self, int bodies);
void simulation_step(body_t *body, int self, int bodies
#ifdef VISUALIZE_OUTPUT
    , int visualize_round
#endif
    );

void send_body(body_t *body, int dest, int tag);
void recv_body(body_t *body, int source, int tag);

int main(int argc, char *argv[]) {
    MPI_Init(&argc, &argv);

    int rank, size;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);

    if(argc != 3) {
        if(rank == 0) printf(USAGE, argv[0]);
        return 1;
    }
    else {
        int rounds = strtol(argv[1], NULL, 0);
        if(rank == 0) {
            body_t *data = parse_input(argv[2], size);
            setup_simulation(data, rounds, 0, size);  // rank == 0
            free(data);
        }
        else {
            setup_simulation(NULL, rounds, rank, size);
        }
    }

    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 self, int bodies) {
    body_t body;

#ifdef RING_COMM_INIT  // ring init
    if(self == 0) {
        for(int i = 1; i < bodies; i ++) {
            send_body(&data[i], self + 1, TAG_INIT);
        }
        body = data[0];
    }
    else {
        recv_body(&body, self - 1, TAG_INIT);
        for(int i = self + 1; i < bodies; i ++) {
            body_t temp_body;
            recv_body(&temp_body, self - 1, TAG_INIT);
            send_body(&temp_body, self + 1, TAG_INIT);
        }
    }
#else  // broadcast init
    if(self == 0) {
        for(int i = 1; i < bodies; i ++) {
            send_body(&data[i], i, TAG_INIT);
        }
        body = data[0];
    }
    else {
        recv_body(&body, 0, TAG_INIT);
    }
#endif

    DEBUG("Start: Node %d has body (%f,%f) (%f,%f) %f\n", self,
        body.pos.x, body.pos.y,
        body.velocity.x, body.velocity.y, body.mass);

    perform_simulation(&body, rounds, self, bodies);

    DEBUG("End: Node %d has body (%f,%f) (%f,%f) %f\n", self,
        body.pos.x, body.pos.y,
        body.velocity.x, body.velocity.y, body.mass);

    DEBUG("ans %d %.9f %.9f %f\n", self, body.pos.x, body.pos.y, body.mass);

#ifdef RING_COMM_TERMINATE  // ring communication
    if(self == 0) {
        printf("%.10f %.10f %f\n", body.pos.x, body.pos.y, body.mass);
        for(int i = 1; i < bodies; i ++) {
            body_t temp_body;
            recv_body(&temp_body, 1, TAG_TERMINATE);
            printf("%.10f %.10f %f\n", temp_body.pos.x, temp_body.pos.y, temp_body.mass);
        }
    }
    else {
        send_body(&body, self - 1, TAG_TERMINATE);
        for(int i = self + 1; i < bodies; i ++) {
            body_t temp_body;
            recv_body(&temp_body, self + 1, TAG_TERMINATE);
            send_body(&temp_body, self - 1, TAG_TERMINATE);
        }
    }
#else  // broadcast communication
    if(self == 0) {
        printf("%.10f %.10f %f\n", body.pos.x, body.pos.y, body.mass);
        for(int i = 1; i < bodies; i ++) {
            body_t temp_body;
            recv_body(&temp_body, i, TAG_TERMINATE);
            printf("%.10f %.10f %f\n", temp_body.pos.x, temp_body.pos.y, temp_body.mass);
        }
    }
    else {
        send_body(&body, 0, TAG_TERMINATE);
    }
#endif
}

void perform_simulation(body_t *body, int rounds, int self, int bodies) {
#ifdef VISUALIZE_OUTPUT
    if(self == 0) printf("%d\n", bodies);
#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, self, bodies, self == 0);
        }
        else {
            simulation_step(body, self, bodies, 0);
        }
#else
        simulation_step(body, self, bodies);
#endif

        DEBUG("Middle: Node %d has body (%f,%f) (%f,%f) %f\n", self,
            body->pos.x, body->pos.y,
            body->velocity.x, body->velocity.y, body->mass);
    }
}

void simulation_step(body_t *body, int self, int bodies
#ifdef VISUALIZE_OUTPUT
    , int visualize_round
#endif
) {
     
    body_t new_body = *body;

    SCALAR force_x = 0;
    SCALAR force_y = 0;
    const SCALAR TIME_DELTA = 1.0;
    const SCALAR CONSTANT_G = 6.6738480e-11;  // from empirical physics

#ifdef RING_COMM_SIMULATE  // ring communication
#ifdef VISUALIZE_OUTPUT
    if(visualize_round) printf("%f %f %f\n", body->pos.x, body->pos.y, body->mass);
#endif
    int prev = (self + bodies - 1) % bodies;
    int next = (self + 1) % bodies;
    int p;
    for(p = 0; p < bodies; p ++) {
        if(p == 0) {
            send_body(body, next, TAG_SIMULATE);
            continue;
        }

        body_t other;
        recv_body(&other, prev, TAG_SIMULATE);
#ifdef VISUALIZE_OUTPUT
        if(visualize_round) printf("%f %f %f\n", other.pos.x, other.pos.y, other.mass);
#endif
        if(p + 1 < bodies) send_body(&other, next, TAG_SIMULATE);
        DEBUG("found at %d: body (%f,%f) (%f,%f) %f\n", self,
            body->pos.x, body->pos.y,
            body->velocity.x, body->velocity.y, body->mass);

        vector_t pos_vector;
        pos_vector.x = other.pos.x - body->pos.x;
        pos_vector.y = other.pos.y - body->pos.y;

        SCALAR length = HYPOT(pos_vector.x, pos_vector.y);
        SCALAR force = (CONSTANT_G * body->mass * other.mass)
            / (length * length);
        force_x += force * (pos_vector.x / length);
        force_y += force * (pos_vector.y / length);
    }
#else  // broadcast communication
    int source, dest;
    for(source = 0; source < bodies; source ++) {
        if(source == self) {
            for(dest = 0; dest < bodies; dest ++) {
                if(dest == self) continue;
                send_body(&new_body, dest, TAG_SIMULATE);
            }
            continue;
        }

        body_t other;
        recv_body(&other, source, TAG_SIMULATE);

        vector_t pos_vector;
        pos_vector.x = other.pos.x - body->pos.x;
        pos_vector.y = other.pos.y - body->pos.y;

        SCALAR length = HYPOT(pos_vector.x, pos_vector.y);
        SCALAR force = (CONSTANT_G * body->mass * other.mass)
            / (length * length);
        force_x += force * (pos_vector.x / length);
        force_y += force * (pos_vector.y / length);
    }
#endif

    new_body.velocity.x = body->velocity.x + force_x * TIME_DELTA / body->mass;
    new_body.velocity.y = body->velocity.y + force_y * TIME_DELTA / body->mass;
    new_body.pos.x = body->pos.x + new_body.velocity.x * TIME_DELTA;
    new_body.pos.y = body->pos.y + new_body.velocity.y * TIME_DELTA;

    *body = new_body;
}

void send_body(body_t *body, int dest, int tag) {
#ifndef SERIALIZE_VELOCITIES
    SCALAR array[] = { body->pos.x, body->pos.y, body->mass };
#else
    SCALAR array[] = {
        body->pos.x, body->pos.y,
        body->velocity.x, body->velocity.y,
        body->mass
    };
#endif
    MPI_Send(array, sizeof(array) / sizeof(*array), SCALAR_MPI_TYPE,
        dest, tag, MPI_COMM_WORLD);
}

void recv_body(body_t *body, int source, int tag) {
#ifndef SERIALIZE_VELOCITIES
    SCALAR array[3];
#else
    SCALAR array[5];
#endif
    MPI_Status status;
    
    MPI_Recv(array, sizeof(array) / sizeof(*array), SCALAR_MPI_TYPE,
        source, tag, MPI_COMM_WORLD, &status);
    
    int i = 0;
    body->pos.x = array[i++];
    body->pos.y = array[i++];
#ifndef SERIALIZE_VELOCITIES
    body->velocity.x = 0.0;
    body->velocity.y = 0.0;
#else
    body->velocity.x = array[i++];
    body->velocity.y = array[i++];
#endif
    body->mass = array[i++];
}
Page generated on Tue Oct 24 00:37:52 2017