kmazouzi / FormationMeso

Download as
Browse Code ยป

Commit 1dea95fb9241ca7e1e9b578aca2b2d0aced9efdf

Authored by kmazouzi
1 parent 3604dfae01
Exists in master

add upper border to 100

Showing 1 changed file with 101 additions and 103 deletions Inline Diff

lab3/mpi_heat2D.c View file @ 1dea95f
/**************************************************************************** 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
printf ( "\n" ); 127
printf ( "HEATED_PLATE\n" ); 128
printf ( " Parallel MPI version using %d processors \n",size ); 129
printf ( " A program to solve for the steady state temperature distribution\n" ); 130
printf ( " over a rectangular plate.\n" ); 131
printf ( "\n" ); 132
printf ( " Spatial grid of %d by %d points.\n", N, N ); 133
printf ( " Each processor will use grid of %d +2 by %d points.\n", M-2, N ); 134
printf ( " The iteration will end until tolerance <= %f\n\n",EPSILON); 135
136
} 137
138
/* Initialize grid and create input file 139
* each process initialize its part 140
* */ 141
142 126
inidat(rank,size,M,N,u,unew); 143 127 printf ( "\n" );
128 printf ( "HEATED_PLATE\n" );
129 printf ( " Parallel MPI version using %d processors \n",size );
130 printf ( " A program to solve for the steady state temperature distribution\n" );
131 printf ( " over a rectangular plate.\n" );
132 printf ( "\n" );
133 printf ( " Spatial grid of %d by %d points.\n", N, N );
134 printf ( " Each processor will use grid of %d by %d points.\n", M, N );
135 printf ( " The iteration will end until tolerance <= %f\n\n",EPSILON);
144 136
prtdat(rank,size,M,N,u, "initial.dat"); 145 137 }
146 138
139 /* Initialize grid and create input file
140 * each process initialize its part
141 * */
147 142
/* 148 143 inidat(rank,size,M,N,u,unew);
* iterate until the new solution unew differs from the old solution u 149 144
145 prtdat(rank,size,M,N,u, "initial.dat");
146
147
148 float diff, globaldiff=1.0;
149 int iter=0;
150
151 /*
152 * iterate (JACOBI ITERATION) until the new solution unew differs from the old solution u
* by no more than EPSILON. 150 153 * by no more than EPSILON.
* */ 151 154 **/
152 155
float diff=1.0; 153 156 while(globaldiff> EPSILON) {
int iter=0; 154
155 157
while(diff> EPSILON) { 156 158 diff= update(rank,size,M,N, u, unew);
157 159
diff= update(rank,size,M,N, u, unew); 158 160 /**
161 * COMPUTE GLOBAL CONVERGENCE
162 * */
159 163
164 MPI_Allreduce(&diff, &globaldiff , 1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);
165
166
if(rank==0) 160 167 if(rank==0)
if(iter%ITER_PRINT==0) 161 168 if(iter%ITER_PRINT==0)
printf("Processor #%d, iteration %d, epsilon = %f\n ", rank,iter,diff); 162 169 printf("Processor #%d, iteration %d, epsilon = %f\n ", rank,iter,globaldiff);
iter++; 163 170 iter++;
} 164 171 }
165 172
prtdat(rank,size,M,N, u, "final.dat"); 166 173 prtdat(rank,size,M,N, u, "final.dat");
free(u); 167 174 free(u);
free(unew); 168 175 free(unew);
MPI_Finalize(); 169 176 MPI_Finalize();
} 170 177 }
171 178
172 179
173 180
/**************************************************************************** 174 181 /****************************************************************************
* subroutine update 175 182 * subroutine update
****************************************************************************/ 176 183 ****************************************************************************/
float update(int rank, int size, int nx,int ny, float *u, float *unew){ 177 184 float update(int rank, int size, int nx,int ny, float *u, float *unew){
int ix, iy; 178 185 int ix, iy;
float diff=0.0; 179 186 float diff=0.0;
float globaldiff; 180
MPI_Status status; 181 187 MPI_Status status;
182 188
183 189
/* 184 190 /*
* EXCHANGE GHOST CELL 185 191 * EXCHANGE GHOST CELL
*/ 186 192 */
if (rank > 0 && rank< size-1) 187 193 if (rank > 0 && rank< size-1)
{ 188 194 {
MPI_Sendrecv(&u[ny*(nx-2)], ny, MPI_FLOAT, rank+1, 0, 189 195 MPI_Sendrecv(&u[ny*(nx-2)], ny, MPI_FLOAT, rank+1, 0,
&u[ny*0], ny, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &status); 190 196 &u[ny*0], ny, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &status);
MPI_Sendrecv(&u[ny*1], ny, MPI_FLOAT, rank-1, 1, 191 197 MPI_Sendrecv(&u[ny*1], ny, MPI_FLOAT, rank-1, 1,
&u[ny*(nx-1)], ny, MPI_FLOAT, rank+1, 1, MPI_COMM_WORLD, &status); 192 198 &u[ny*(nx-1)], ny, MPI_FLOAT, rank+1, 1, MPI_COMM_WORLD, &status);
} 193 199 }
194 200
else if (rank == 0 && rank< size-1) 195 201 else if (rank == 0)
MPI_Sendrecv(&u[ny*(nx-2)], ny, MPI_FLOAT, rank+1, 0, 196 202 MPI_Sendrecv(&u[ny*(nx-2)], ny, MPI_FLOAT, rank+1, 0,
&u[ny*(nx-1)], ny, MPI_FLOAT, rank+1, 1, MPI_COMM_WORLD, &status); 197 203 &u[ny*(nx-1)], ny, MPI_FLOAT, rank+1, 1, MPI_COMM_WORLD, &status);
else if ( rank> 0 && rank == size-1) 198 204 else if (rank == size-1)
MPI_Sendrecv(&u[ny*1], ny, MPI_FLOAT, rank-1, 1, 199 205 MPI_Sendrecv(&u[ny*1], ny, MPI_FLOAT, rank-1, 1,
&u[ny*0], ny, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &status); 200 206 &u[ny*0], ny, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &status);
201 207
202 208
203 209
/** 204 210 /**
* PERFORM LOCAL COMPUTATION 205 211 * PERFORM LOCAL COMPUTATION
* */ 206 212 * */
207 213
for (ix = 1; ix < nx-1; ix++) { 208 214 for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) { 209 215 for (iy = 1; iy < ny-1; iy++) {
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 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
; 211 217 ;
if (diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] )) 212 218 if (diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] ))
{ 213 219 {
diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] ); 214 220 diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] );
} 215 221 }
} 216 222 }
217 223
} 218 224 }
219 225
220 226
/** 221
* COMPUTE GLOBAL CONVERGENCE 222
* 223
* */ 224
225
MPI_Allreduce(&diff, &globaldiff , 1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD); 226
227
228
/** 229
* COPY OLD DATA 230
* */ 231
232
233
for (ix = 1; ix < nx-1; ix++) { 234 227 for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) { 235 228 for (iy = 1; iy < ny-1; iy++) {
u[ix*ny+iy] = unew[ix*ny+iy]; 236 229 u[ix*ny+iy] = unew[ix*ny+iy];
} 237 230 }
} 238 231 }
239 232
240 233
return globaldiff; 241 234 return diff;
} 242 235 }
243 236
/***************************************************************************** 244 237 /*****************************************************************************
* Initialize Data 245 238 * Initialize Data
*****************************************************************************/ 246 239 *****************************************************************************/
247 240
void inidat(int rank, int size,int nx, int ny, float *u, float *unew) 248 241 void inidat(int rank, int size,int nx, int ny, float *u, float *unew)
{ 249 242 {
int ix, iy; 250 243 int ix, iy;
251 244
/* 252 245 /*
*Set boundary data and interrior values 253 246 *Set boundary data and interrior values
* */ 254 247 * */
255 248
256 249
// interior points 257 250 // interior points
for (ix = 1; ix < nx-1; ix++) 258 251 for (ix = 1; ix < nx-1; ix++)
for (iy = 1; iy < ny-1; iy++) { 259 252 for (iy = 1; iy < ny-1; iy++) {
u[ix*ny+iy]=0.0; 260 253 u[ix*ny+iy]=0.0;
} 261 254 }
262 255
//boundary left 263 256 //boundary left
for (ix = 1; ix < nx-1; ix++){ 264 257 for (ix = 0; ix < nx; ix++){
u[ix*ny]=100.0; 265 258 u[ix*ny]=100.0;
266 259
} 267 260 }
268 261
//boundary right 269 262 //boundary right
for (ix = 1; ix < nx-1; ix++){ 270 263 for (ix = 0; ix < nx; ix++){
u[ix*ny+ (ny-1)]=100.0; 271 264 u[ix*ny+ (ny-1)]=100.0;
272 265
} 273 266 }
274 267
//boundary down 275 268 //boundary down
for (iy = 0; iy < ny; iy++){ 276 269 for (iy = 1; iy < ny-1; iy++){
277 270
if(rank==size-1) { 278 271 if(rank==size-1) {
u[(nx-1)*(ny)+iy]=100.0; 279 272 u[(nx-1)*(ny)+iy]=100.0;
}else 280 273 u[0]=100.0;
{ 281 274
u[(nx-1)*(ny)+iy]=0.0; 282 275 }else
} 283 276 {
277 u[(nx-1)*(ny)+iy]=0.0;
278 }
} 284 279 }
285 280
//boundary top 286 281 //boundary top
for (iy = 0; iy < ny; iy++){ 287 282 for (iy = 1; iy < ny-1; iy++){
u[iy]=0.0; 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
294 292
/*************************************************************************** 295 293 /***************************************************************************
* Print Data to files 296 294 * Print Data to files
**************************************************************************/ 297 295 **************************************************************************/
298 296
void print2file(int rank, int nx, int ny, float *u,const char *fname) 299 297 void print2file(int rank, int nx, int ny, float *u,const char *fname)
{ 300 298 {
int ix, iy; 301 299 int ix, iy;
FILE *fp; 302 300 FILE *fp;
303 301
char str[255]; 304 302 char str[255];
305 303
sprintf(str, "%d_%s", rank,fname); 306 304 sprintf(str, "%d_%s", rank,fname);