Skip to content

Instantly share code, notes, and snippets.

@michaelchum

michaelchum/grid_512_512.c

Created November 11, 2016 21:30
Show Gist options
  • Save michaelchum/da2870fb0de40f4f50be64975733acad to your computer and use it in GitHub Desktop.
Save michaelchum/da2870fb0de40f4f50be64975733acad to your computer and use it in GitHub Desktop.
Download ZIP
Raw
grid_512_512.c
#include <mpi.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <time.h>
#define G 0.75
#define N 512
#define TAG 0
#define ETA 0.0002
#define RHO 0.5
int main(int argc, char * argv[]) {
// Number of iterations
int iterations = atoi(argv[1]);
// Initialize MPI
MPI_Init(&argc, &argv);
// Get current rank and number of processes
int rank, size;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
// Make a first_comm
int first_color = rank/2;
MPI_Comm first_comm;
MPI_Comm_split(MPI_COMM_WORLD, first_color, rank, &first_comm);
// Make a second_comm
int second_color = (rank-1)/2;
MPI_Comm second_comm;
MPI_Comm_split(MPI_COMM_WORLD, second_color, rank, &second_comm);
int first_comm_rank, first_comm_size;
MPI_Comm_rank(first_comm, &first_comm_rank);
MPI_Comm_size(first_comm, &first_comm_size);
int second_comm_rank, second_comm_size;
MPI_Comm_rank(second_comm, &second_comm_rank);
MPI_Comm_size(second_comm, &second_comm_size);
// Initialize 3 matrices, one for each time state
float u[N][N] = {0};
float u1[N][N] = {0};
float u2[N][N] = {0};
// Initialize state for N/2, N/2 e.g. Simulate drum hit!
u1[N/2][N/2] += 1;
// Assume size if even or 1
int block_size = N * N / size;
int num_rows_each = N / size;
int start = rank * num_rows_each;
int end = num_rows_each + start -1;
MPI_Status status;
FILE * fp;
// Print to file
if(rank==(size/2)){
fp = fopen("output_512.h","w");
fprintf(fp,"float output[%d] = {\n", iterations);
}
for (int i = 0; i < iterations; i++){
if(rank == 0){
MPI_Sendrecv(&u1[end][0], N, MPI_FLOAT, first_comm_rank+1, TAG, &u1[end+1][0], N, MPI_FLOAT, first_comm_rank+1, TAG, first_comm, &status);
} else if (rank == size-1) {
MPI_Sendrecv(&u1[start][0], N, MPI_FLOAT, first_comm_rank-1, TAG, &u1[start-1][0], N, MPI_FLOAT, first_comm_rank-1, TAG, first_comm, &status);
} else if(rank%2 == 0) {
MPI_Sendrecv(&u1[end][0], N, MPI_FLOAT, first_comm_rank+1, TAG, &u1[end+1][0], N, MPI_FLOAT, first_comm_rank+1, TAG, first_comm, &status);
MPI_Sendrecv(&u1[start][0], N, MPI_FLOAT, second_comm_rank-1, TAG, &u1[start-1][0], N, MPI_FLOAT, second_comm_rank-1, TAG, second_comm, &status);
} else {
MPI_Sendrecv(&u1[start][0], N, MPI_FLOAT, first_comm_rank-1, TAG, &u1[start-1][0], N, MPI_FLOAT, first_comm_rank-1, TAG, first_comm, &status);
MPI_Sendrecv(&u1[end][0], N, MPI_FLOAT, second_comm_rank+1, TAG, &u1[end+1][0], N, MPI_FLOAT, second_comm_rank+1, TAG, second_comm, &status);
}
// Wait for all communications to be done across all ranks
MPI_Barrier(MPI_COMM_WORLD);
// CENTER
for (int i = 1; i <= N-2; i++) {
if (i>=start && i<=end){
int k;
for(k = 1; k <= N-2; k++){
u[i][k] = ((RHO) * (u1[i-1][k] + u1[i+1][k] + u1[i][k-1] + u1[i][k+1] - (4 * u1[i][k])) + 2 * u1[i][k] - (1 - ETA) * u2[i][k])/(1 + ETA);
}
}
}
// SIDES
for (int i = 1; i <= N - 2; i++) {
u[0][i] = G * u[1][i];
u[N-1][i] = G * u[N-2][i];
u[i][0] = G * u[i][1];
u[i][N-1] = G * u[i][N-2];
}
// CORNERS
u[0][0] = G * u[1][0];
u[N-1][0] = G * u[N-2][0];
u[0][N-1] = G * u[0][N-2];
u[N-1][N-1] = G * u[N-1][N-2];
if (rank==(size/2)) {
if (i==iterations-1) {
printf("%f\n", u[N/2][N/2]);
// Print to file
fprintf(fp,"%f\n",u[N/2][N/2]);
fprintf(fp,"};\n");
fclose(fp);
}
else {
printf("%f,\n", u[N/2][N/2]);
// Print to file
fprintf(fp,"%f,\n",u[N/2][N/2]);
}
}
// Update the matrices
memcpy(u2, u1, sizeof(u2));
memcpy(u1, u, sizeof(u1));
}
MPI_Comm_free(&first_comm);
MPI_Comm_free(&second_comm);
MPI_Finalize();
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment