/******************************************************************************* Erik Barry Erhardt and Matthew Bohnsack CS 442, Project 3, 4/30/2006 wave.c qsub -I -l nodes=4,walltime=3600 touch wave.c make mpirun -np 16 -machinefile $PBS_NODEFILE ./wave 2 240 4 4 *******************************************************************************/ #include #include #include #define MAX_DIM 800 int rank, rank_x, rank_y, rank_down, rank_up, rank_left, rank_right; int nprocs, nprocsx, nprocsy; int nxg, nyg, nx, ny, n; double length_x, length_y, delta_x, delta_y; double PI, h, t0, tfinal, t, time_per_step; int junk, error_flag=0, sw_method; // this is our data matrix double U[MAX_DIM+2][MAX_DIM+2], Ue[MAX_DIM+2][MAX_DIM+2], lu[MAX_DIM+2][MAX_DIM+2]; double U1[MAX_DIM+2][MAX_DIM+2]; // this is our data at time n-1 /******************************************************************************/ int main(int argc, char *argv[]){ int i, j, it; double time1, time2, timetotal, timesofar, timetogo; // time the run MPI_Init(&argc,&argv); // initialize MPI MPI_Comm_rank(MPI_COMM_WORLD, &rank); // my rank number MPI_Comm_size(MPI_COMM_WORLD, &nprocs); // number of processes time1 = MPI_Wtime(); // first time t1 junk = init(argc, &*argv); /*******/ // initialize values if(sw_method==2){ for (it=0; it<=n+1; it++){ if(it % ((n-1)/10) == 0) { // output grid to a file double x_out[nx]; double y_out[ny]; double u_out[nx][ny]; for (i=0; i= pow(delta_x,2)/8){if(rank==0){printf("ERROR: h >= pow(delta_x,2)/8 (%e >= %e ), // EXIT\n", h, pow(delta_x,2)/8);} error_flag++;} if(error_flag>0){if(rank==0) printf("Errors = %d\n EXITING\n", error_flag); exit(0);} junk = init_local_ranks(); junk = init_block_data(); return 0; } int init_local_ranks(){ // if(rank==0) printf("nprocs x y %d %d %d\n", nprocs, nprocsx, nprocsy); // x, y rank location of current process rank_x = rank % nprocsx; // rank in x direction (mod) rank_y = rank / nprocsx; // rank in y direction (int division) /* Ranks (blocks) go in this order 6 7 8 3 4 5 0 1 2 */ // ranks of neighboring processors rank_down = (rank_y-1)*nprocsx + rank_x ; if(rank_y == 0 ) rank_down = -1; rank_up = (rank_y+1)*nprocsx + rank_x ; if(rank_y == nprocsy-1) rank_up = -1; rank_left = rank_y *nprocsx + rank_x-1; if(rank_x == 0 ) rank_left = -1; rank_right = rank_y *nprocsx + rank_x+1; if(rank_x == nprocsx-1) rank_right = -1; //printf("Ranks: me, x, y, up, down, left, right: %d %d %d %d %d %d %d \n", rank, rank_x, // rank_y, rank_up, rank_down, rank_left, rank_right); return 0; } int init_block_data(){ int i, j; double b, initx[nx+2], inity[nx+2], pi2; /* Data goes in this order per block (border is ghost cells): 0 1 2 3 4 ... nx+1 1 2 ... ... ny+1 */ for (i=0; i0) { // simple scheme term1 = U[i-1][j] - 2*U[i][j] + U[i+1][j]; term2 = U[i][j-1] - 2*U[i][j] + U[i][j+1]; term3 = 2*U[i][j]-U1[i][j]; lu[i][j] = pow(a,2)*pow(h,2)*(term1/pow(delta_x,2)+term2/pow(delta_y,2))+term3; } else { // initial time step term1 = U[i-1][j] - U[i+1][j]; term2 = U[i][j-1] - U[i][j+1]; term3 = U[i][j]; lu[i][j] = a*h*(term1/delta_x+term2/delta_y)+term3; } } } } return 0; } /******************************************************************************/ /* ghost_cell_update *****************************************/ int ghost_cell_update(){ int i, j; int ierr, ierr1, ierr2, ierr3, tag=0; MPI_Status status_rreq[4], status_sreq[4]; MPI_Request rreq[4], sreq[4]; // receive and send requests (down, up, left, right) // receive and send buffers double bufrd[nx],bufsd[nx],bufru[nx],bufsu[nx],bufrl[ny],bufsl[ny],bufrr[ny],bufsr[ny]; // DOWN ********************** if(rank_down==-1){ // down boundary condition; //j=ny; for (i=1; i<=nx; i++){U[i][j]=U[i][j+1];} sreq[0] = MPI_SUCCESS; rreq[0] = MPI_SUCCESS; } else { // MPI j=ny; for (i=1; i= '0' && *digit <='9') { result = (result * 10) + (*digit - '0'); digit++; } //--- Check that there were no non-digits at end. if (*digit != 0) {return -1;} return result; } /******************************************************************************/ /* output *****************************************/ int output(int nx,int ny,double x[],double y[],double u[nx][ny],double time) { double xout[nx],yout[ny],uout[nx][ny]; double uout2[ny][nx]; double tmp; int ierr,amode,rank,nproc,p,i,j; MPI_File fh; static int firstcall=1; char *fname="wave.out"; MPI_Comm_rank(MPI_COMM_WORLD,&rank); MPI_Comm_size(MPI_COMM_WORLD,&nproc); if(rank==0) printf("Write: %s at t=%f\n",fname,time); if (firstcall) { amode = MPI_MODE_WRONLY | MPI_MODE_CREATE; MPI_File_delete(fname,MPI_INFO_NULL); firstcall=0; }else{ amode = MPI_MODE_WRONLY | MPI_MODE_APPEND; } MPI_File_open(MPI_COMM_WORLD,fname,amode,MPI_INFO_NULL,&fh); if (rank==0) { double tmp=nproc; MPI_File_write(fh,&tmp,1,MPI_DOUBLE,MPI_STATUSES_IGNORE); tmp=nx; MPI_File_write(fh,&tmp,1,MPI_DOUBLE,MPI_STATUSES_IGNORE); tmp=ny; MPI_File_write(fh,&tmp,1,MPI_DOUBLE,MPI_STATUSES_IGNORE); } for (p=0; p