kmazouzi / FormationMeso

Compare View

switch
from
...
to
 
Commits (2)

Diff

Showing 2 changed files Inline Diff

lab3/concat.sh View file @ 6c0f7fa
File was created 1 #!/bin/bash
2 rm -rf initial.dat final.dat
3 cat *_initial.dat >> initial.dat
4 cat *_final.dat >> final.dat
lab3/mpi_heat2D.c View file @ 6c0f7fa
/**************************************************************************** 1 1 /****************************************************************************
* DESCRIPTION: 2 2 * DESCRIPTION:
* Serial HEAT2D Example - C Version 3 3 * Serial HEAT2D Example - C Version
* This example is based on a simplified 4 4 * This example is based on a simplified
* two-dimensional heat equation domain decomposition. The initial 5 5 * two-dimensional heat equation domain decomposition. The initial
* temperature is computed to be high in the middle of the domain and 6 6 * temperature is computed to be high in the middle of the domain and
* zero at the boundaries. The boundaries are held at zero throughout 7 7 * zero at the boundaries. The boundaries are held at zero throughout
* the simulation. During the time-stepping, an array containing two 8 8 * the simulation. During the time-stepping, an array containing two
* domains is used; these domains alternate between old data and new data. 9 9 * domains is used; these domains alternate between old data and new data.
* 10 10 *
* The physical region, and the boundary conditions, are suggested 11 11 * The physical region, and the boundary conditions, are suggested
by this diagram; 12 12 by this diagram;
13 13
u = 0 14 14 u = 0
+------------------+ 15 15 +------------------+
| | 16 16 | |
u = 100 | | u = 100 17 17 u = 100 | | u = 100
| | 18 18 | |
| | 19 19 | |
| | 20 20 | |
| | 21 21 | |
+------------------+ 22 22 +------------------+
u = 100 23 23 u = 100
24 24
Interrior point : 25 25 Interrior point :
u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] ) 26 26 u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] )
27 27
28 28
PARALLEL MPI VERSION : 29 29 PARALLEL MPI VERSION :
30 30
+-------------------+ 31 31 +-------------------+
| | P0 m=(n-2)/P +2 32 32 | | P0 m=(n-2)/P +2
+-------------------+ 33 33 +-------------------+
| | P1 34 34 | | P1
+-------------------+ 35 35 +-------------------+
n | | .. 36 36 n | | ..
+-------------------+ 37 37 +-------------------+
| | Pq 38 38 | | Pq
+-------------------+ 39 39 +-------------------+
40 40
<-------- n --------> 41 41 <-------- n -------->
<-------n-2 ------> 42 42 <-------n-2 ------>
43 43
44 44
****************************************************************************/ 45 45 ****************************************************************************/
#include <stdio.h> 46 46 #include <stdio.h>
#include <stdlib.h> 47 47 #include <stdlib.h>
#include <math.h> 48 48 #include <math.h>
#include <mpi/mpi.h> 49 49 #include <mpi/mpi.h>
#define NN 50 50 50 #define NN 50
#define MM 50 51 51 #define MM 50
52 52
#define RING 100 53 53 #define RING 100
#define ITER_PRINT 100 54 54 #define ITER_PRINT 100
#define PRINT_DATA 1 55 55 #define PRINT_DATA 1
#define MAX_ITER 1000 56 56 #define MAX_ITER 1000
#define _EPSILON 0.01 57 57 #define _EPSILON 0.01
58 58
59 59
float update(int rank,int size, int nx,int ny, float *u, float *unew); 60 60 float update(int rank,int size, int nx,int ny, float *u, float *unew);
void inidat(int rank, int size, int nx, int ny, float *u, float *unew); 61 61 void inidat(int rank, int size, int nx, int ny, float *u, float *unew);
void prtdat(int rank, int size, int nx, int ny, float *u,const char *fnam); 62 62 void prtdat(int rank, int size, int nx, int ny, float *u,const char *fnam);
63 63
64 64
65 65
66 66
int main(int argc, char *argv[]) 67 67 int main(int argc, char *argv[])
{ 68 68 {
int N=NN,M=MM; 69 69 int N=NN,M=MM;
70 70
int rank,size; 71 71 int rank,size;
72 72
float EPSILON=_EPSILON; 73 73 float EPSILON=_EPSILON;
74 74
75 75
/* INITIALIZE MPI */ 76 76 /* INITIALIZE MPI */
MPI_Init(&argc, &argv); 77 77 MPI_Init(&argc, &argv);
78 78
/* GET THE PROCESSOR ID AND NUMBER OF PROCESSORS */ 79 79 /* GET THE PROCESSOR ID AND NUMBER OF PROCESSORS */
MPI_Comm_rank(MPI_COMM_WORLD, &rank); 80 80 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size); 81 81 MPI_Comm_size(MPI_COMM_WORLD, &size);
82 82
//Only Rank 0 read application parameters 83 83 //Only Rank 0 read application parameters
if(rank==0) { 84 84 if(rank==0) {
85 85
if(argc !=3) 86 86 if(argc !=3)
{ 87 87 {
fprintf(stderr,"usage %s N EPSILON\n ", argv[0]); 88 88 fprintf(stderr,"usage %s N EPSILON\n ", argv[0]);
fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n"); 89 89 fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n");
fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n"); 90 90 fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n");
return -1; 91 91 return -1;
} 92 92 }
93 93
N = M = atoi(argv[1]); 94 94 N = M = atoi(argv[1]);
EPSILON = atof(argv[2]); 95 95 EPSILON = atof(argv[2]);
96 96
if(N % size!=0) 97 97 if(N % size!=0)
{ 98 98 {
fprintf(stderr,"Grid Size MUST be divisible by the number of processors !"); 99 99 fprintf(stderr,"Grid Size MUST be divisible by the number of processors !");
return -1; 100 100 return -1;
} 101 101 }
102 102
} 103 103 }
104 104
//Wait for rank 0 , all process start here 105 105 //Wait for rank 0 , all process start here
MPI_Barrier(MPI_COMM_WORLD); 106 106 MPI_Barrier(MPI_COMM_WORLD);
107 107
//Exchange N 108 108 //Exchange N
MPI_Bcast(&N , 1, MPI_FLOAT, 0 , MPI_COMM_WORLD); 109 109 MPI_Bcast(&N , 1, MPI_FLOAT, 0 , MPI_COMM_WORLD);
//Exchange EPSILON 110 110 //Exchange EPSILON
MPI_Bcast(&EPSILON , 1, MPI_FLOAT, 0 , MPI_COMM_WORLD); 111 111 MPI_Bcast(&EPSILON , 1, MPI_FLOAT, 0 , MPI_COMM_WORLD);
112 112
//local size 113 113 //local size
M = (N-2) / size + 2 ; 114 114 M = (N-2) / size + 2 ;
115 115
float *u = (float *)malloc(N * M * sizeof(float)); 116 116 float *u = (float *)malloc(N * M * sizeof(float));
float *unew = (float *)malloc(N * M * sizeof(float)); 117 117 float *unew = (float *)malloc(N * M * sizeof(float));
118 118
if(u==0 || unew ==0) 119 119 if(u==0 || unew ==0)
{ 120 120 {
perror("Can't allocated data\n"); 121 121 perror("Can't allocated data\n");
return -1; 122 122 return -1;
} 123 123 }
124 124
if(rank==0) { 125 125 if(rank==0) {
126 126
printf ( "\n" ); 127 127 printf ( "\n" );
printf ( "HEATED_PLATE\n" ); 128 128 printf ( "HEATED_PLATE\n" );
printf ( " Parallel MPI version using %d processors \n",size ); 129 129 printf ( " Parallel MPI version using %d processors \n",size );
printf ( " A program to solve for the steady state temperature distribution\n" ); 130 130 printf ( " A program to solve for the steady state temperature distribution\n" );
printf ( " over a rectangular plate.\n" ); 131 131 printf ( " over a rectangular plate.\n" );
printf ( "\n" ); 132 132 printf ( "\n" );
printf ( " Spatial grid of %d by %d points.\n", N, N ); 133 133 printf ( " Spatial grid of %d by %d points.\n", N, N );
printf ( " Each processor will use grid of %d by %d points.\n", M, N ); 134 134 printf ( " Each processor will use grid of %d by %d points.\n", M, N );
printf ( " The iteration will end until tolerance <= %f\n\n",EPSILON); 135 135 printf ( " The iteration will end until tolerance <= %f\n\n",EPSILON);
136 136
} 137 137 }
138 138
/* Initialize grid and create input file 139 139 /* Initialize grid and create input file
* each process initialize its part 140 140 * each process initialize its part
* */ 141 141 * */
142 142
inidat(rank,size,M,N,u,unew); 143 143 inidat(rank,size,M,N,u,unew);
144 144
prtdat(rank,size,M,N,u, "initial.dat"); 145 145 prtdat(rank,size,M,N,u, "initial.dat");
146 146
147 147
float diff, globaldiff=1.0; 148 148 float diff, globaldiff=1.0;
int iter=0; 149 149 int iter=0;
150
151 /*
152 * iterate (JACOBI ITERATION) until the new solution unew differs from the old solution u
150 153 * by no more than EPSILON.
/* 151 154 **/
155
156 while(globaldiff> EPSILON) {
* iterate (JACOBI ITERATION) until the new solution unew differs from the old solution u 152 157
* by no more than EPSILON. 153 158 diff= update(rank,size,M,N, u, unew);
**/ 154
155 159
while(globaldiff> EPSILON) { 156 160 /**
161 * COMPUTE GLOBAL CONVERGENCE
162 * */
163
164 MPI_Allreduce(&diff, &globaldiff , 1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);
157 165
diff= update(rank,size,M,N, u, unew); 158
159 166
/** 160 167 if(rank==0)
* COMPUTE GLOBAL CONVERGENCE 161 168 if(iter%ITER_PRINT==0)
* */ 162 169 printf("Processor #%d, iteration %d, epsilon = %f\n ", rank,iter,globaldiff);
163 170 iter++;
MPI_Allreduce(&diff, &globaldiff , 1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD); 164 171 }
165 172
166 173 prtdat(rank,size,M,N, u, "final.dat");
if(rank==0) 167 174 free(u);
if(iter%ITER_PRINT==0) 168 175 free(unew);
printf("Processor #%d, iteration %d, epsilon = %f\n ", rank,iter,globaldiff); 169 176 MPI_Finalize();
iter++; 170 177 }
} 171 178
172 179
prtdat(rank,size,M,N, u, "final.dat"); 173 180
free(u); 174 181 /****************************************************************************
free(unew); 175 182 * subroutine update
MPI_Finalize(); 176 183 ****************************************************************************/
} 177 184 float update(int rank, int size, int nx,int ny, float *u, float *unew){
178 185 int ix, iy;
179 186 float diff=0.0;
180
/**************************************************************************** 181 187 MPI_Status status;
* subroutine update 182 188
****************************************************************************/ 183 189
float update(int rank, int size, int nx,int ny, float *u, float *unew){ 184 190 /*
int ix, iy; 185 191 * EXCHANGE GHOST CELL
float diff=0.0; 186 192 */
MPI_Status status; 187 193 if (rank > 0 && rank< size-1)
188 194 {
189 195 MPI_Sendrecv(&u[ny*(nx-2)], ny, MPI_FLOAT, rank+1, 0,
/* 190 196 &u[ny*0], ny, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &status);
* EXCHANGE GHOST CELL 191 197 MPI_Sendrecv(&u[ny*1], ny, MPI_FLOAT, rank-1, 1,
*/ 192 198 &u[ny*(nx-1)], ny, MPI_FLOAT, rank+1, 1, MPI_COMM_WORLD, &status);
if (rank > 0 && rank< size-1) 193 199 }
{ 194 200
MPI_Sendrecv(&u[ny*(nx-2)], ny, MPI_FLOAT, rank+1, 0, 195 201 else if (rank == 0)
&u[ny*0], ny, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &status); 196 202 MPI_Sendrecv(&u[ny*(nx-2)], ny, MPI_FLOAT, rank+1, 0,
MPI_Sendrecv(&u[ny*1], ny, MPI_FLOAT, rank-1, 1, 197 203 &u[ny*(nx-1)], ny, MPI_FLOAT, rank+1, 1, MPI_COMM_WORLD, &status);
&u[ny*(nx-1)], ny, MPI_FLOAT, rank+1, 1, MPI_COMM_WORLD, &status); 198 204 else if (rank == size-1)
} 199 205 MPI_Sendrecv(&u[ny*1], ny, MPI_FLOAT, rank-1, 1,
200 206 &u[ny*0], ny, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &status);
else if (rank == 0) 201 207
MPI_Sendrecv(&u[ny*(nx-2)], ny, MPI_FLOAT, rank+1, 0, 202 208
&u[ny*(nx-1)], ny, MPI_FLOAT, rank+1, 1, MPI_COMM_WORLD, &status); 203 209
else if (rank == size-1) 204 210 /**
MPI_Sendrecv(&u[ny*1], ny, MPI_FLOAT, rank-1, 1, 205 211 * PERFORM LOCAL COMPUTATION
&u[ny*0], ny, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &status); 206 212 * */
207 213
208 214 for (ix = 1; ix < nx-1; ix++) {
209 215 for (iy = 1; iy < ny-1; iy++) {
/** 210 216 unew[ix*ny+iy] = (u[(ix+1)*ny+iy] + u[(ix-1)*ny+iy] + u[ix*ny+iy+1] + u[ix*ny+iy-1] )/4.0
* PERFORM LOCAL COMPUTATION 211 217 ;
* */ 212 218 if (diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] ))
213 219 {
for (ix = 1; ix < nx-1; ix++) { 214 220 diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] );
for (iy = 1; iy < ny-1; iy++) { 215 221 }
unew[ix*ny+iy] = (u[(ix+1)*ny+iy] + u[(ix-1)*ny+iy] + u[ix*ny+iy+1] + u[ix*ny+iy-1] )/4.0 216 222 }
; 217 223
if (diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] )) 218 224 }
{ 219 225
diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] ); 220 226
} 221
} 222
223
} 224
225
226
for (ix = 1; ix < nx-1; ix++) { 227
for (iy = 1; iy < ny-1; iy++) { 228
u[ix*ny+iy] = unew[ix*ny+iy]; 229
} 230
} 231
232
233
return diff; 234 227 for (ix = 1; ix < nx-1; ix++) {
} 235 228 for (iy = 1; iy < ny-1; iy++) {
236 229 u[ix*ny+iy] = unew[ix*ny+iy];
/***************************************************************************** 237 230 }
* Initialize Data 238 231 }
*****************************************************************************/ 239 232
240 233
void inidat(int rank, int size,int nx, int ny, float *u, float *unew) 241 234 return diff;
{ 242 235 }
int ix, iy; 243 236
244 237 /*****************************************************************************
/* 245 238 * Initialize Data
*Set boundary data and interrior values 246 239 *****************************************************************************/
* */ 247 240
248 241 void inidat(int rank, int size,int nx, int ny, float *u, float *unew)
249 242 {
// interior points 250 243 int ix, iy;
for (ix = 1; ix < nx-1; ix++) 251 244
for (iy = 1; iy < ny-1; iy++) { 252 245 /*
u[ix*ny+iy]=0.0; 253 246 *Set boundary data and interrior values
} 254 247 * */
255 248
//boundary left 256 249
for (ix = 0; ix < nx; ix++){ 257 250 // interior points
u[ix*ny]=100.0; 258 251 for (ix = 1; ix < nx-1; ix++)
259 252 for (iy = 1; iy < ny-1; iy++) {
} 260 253 u[ix*ny+iy]=0.0;
261 254 }
//boundary right 262 255
for (ix = 0; ix < nx; ix++){ 263 256 //boundary left
u[ix*ny+ (ny-1)]=100.0; 264 257 for (ix = 0; ix < nx; ix++){
265 258 u[ix*ny]=100.0;
} 266 259
267 260 }
//boundary down 268 261
for (iy = 1; iy < ny-1; iy++){ 269 262 //boundary right
270 263 for (ix = 0; ix < nx; ix++){
if(rank==size-1) { 271 264 u[ix*ny+ (ny-1)]=100.0;
u[(nx-1)*(ny)+iy]=100.0; 272 265
u[0]=100.0; 273 266 }
274 267
}else 275 268 //boundary down
{ 276 269 for (iy = 1; iy < ny-1; iy++){
u[(nx-1)*(ny)+iy]=0.0; 277 270
} 278 271 if(rank==size-1) {
} 279 272 u[(nx-1)*(ny)+iy]=100.0;
280 273 u[0]=100.0;
//boundary top 281 274
for (iy = 1; iy < ny-1; iy++){ 282 275 }else
if(rank==0) 283 276 {
277 u[(nx-1)*(ny)+iy]=0.0;
278 }
u[iy]=100.0; 284 279 }
else 285 280
u[iy]=0.0; 286 281 //boundary top
287 282 for (iy = 1; iy < ny-1; iy++){
} 288 283 if(rank==0)
284 u[iy]=100.0;
285 else
286 u[iy]=0.0;
289 287
} 290 288 }
291 289
292 290 }
/*************************************************************************** 293 291
* Print Data to files 294 292
**************************************************************************/ 295 293 /***************************************************************************
296 294 * Print Data to files
void print2file(int rank, int nx, int ny, float *u,const char *fname) 297 295 **************************************************************************/
{ 298 296
int ix, iy; 299 297 void print2file(int rank, int nx, int ny, float *u,const char *fname)
FILE *fp; 300 298 {
301 299 int ix, iy;
char str[255]; 302 300 FILE *fp;
303 301
sprintf(str, "%d_%s", rank,fname); 304 302 char str[255];
305 303
fp = fopen(str, "w"); 306 304 sprintf(str, "%d_%s", rank,fname);
307 305
for (ix = 0 ; ix < nx; ix++) { 308 306 fp = fopen(str, "w");
for (iy =0; iy < ny; iy++) { 309 307
310 308 for (ix = 0 ; ix < nx; ix++) {
fprintf(fp, "%6.2f ", u[ix*ny+iy]); 311 309 for (iy =0; iy < ny; iy++) {
} 312 310
fputc ( '\n', fp); 313 311 fprintf(fp, "%6.2f ", u[ix*ny+iy]);
} 314 312 }
315 313 fputc ( '\n', fp);
fclose(fp); 316 314 }
317 315
} 318 316 fclose(fp);
319 317
320 318 }
321 319
void prtdat(int rank, int size,int nx, int ny, float *u,const char *fname) 322 320
{ 323 321
324 322 void prtdat(int rank, int size,int nx, int ny, float *u,const char *fname)
325 323 {
if(ITER_PRINT==0)return; 326 324
327 325
328 326 if(ITER_PRINT==0)return;
/** 329 327