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);
}
Page generated on Tue Oct 24 00:36:08 2017