kmazouzi / FormationMeso

Download as
Browse Code ยป

Commit d16e53ac39a72d169d9b7c204be121caa81009aa

Authored by kmazouzi
1 parent ead55f0407
Exists in master

add program arguments

Showing 2 changed files with 107 additions and 90 deletions Inline Diff

lab3/omp_heat2D.c View file @ d16e53a
/**************************************************************************** 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 +------------------+
u = 100 20 20 u = 100
21 21
Interrior point : 22 22 Interrior point :
u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] ) 23 23 u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] )
24 24
****************************************************************************/ 25 25 ****************************************************************************/
#include <stdio.h> 26 26 #include <stdio.h>
#include <stdlib.h> 27 27 #include <stdlib.h>
#include <math.h> 28 28 #include <math.h>
#include <omp.h> 29 29 #include <omp.h>
30 30
#define N 500 31 31 #define NN 50
#define M 500 32 32 #define MM 50
33 33
#define ITER_PRINT 100 34 34 #define ITER_PRINT 100
#define PRINT_DATA 1 35 35 #define PRINT_DATA 1
36 36
#define EPSILON 1e-1 37 37 #define _EPSILON 0.001
38 38
39 39
void update(int nx,int ny, float *u, float *unew, float * diff); 40 40 void update(int nx,int ny, float *u, float *unew, float * diff);
void inidat(int nx, int ny, float *u, float *unew); 41 41 void inidat(int nx, int ny, float *u, float *unew);
void prtdat(int nx, int ny, float *u,const char *fnam); 42 42 void prtdat(int nx, int ny, float *u,const char *fnam);
43 43
44 44
45 45
46 46
int main(int argc, char *argv[]) 47 47 int main(int argc, char *argv[])
{ 48 48 {
49 49
float diff=1.0; 50 50 float diff=1.0;
51 float EPSILON=_EPSILON;
52 int N=NN,M=MM;
51 53
54 if(argc !=3)
55 {
56 fprintf(stderr,"usage %s N EPSILON\n ", argv[0]);
57 fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n");
58 fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n");
59 return -1;
60 }
61
62 N = M = atoi(argv[1]);
63 EPSILON = atof(argv[2]);
64
float *u = (float *)malloc(N * M * sizeof(float)); 52 65 float *u = (float *)malloc(N * M * sizeof(float));
float *unew = (float *)malloc(N * M * sizeof(float)); 53 66 float *unew = (float *)malloc(N * M * sizeof(float));
54 67
if(u==0 || unew ==0) 55 68 if(u==0 || unew ==0)
{ 56 69 {
perror("Can't allocated data\n"); 57 70 perror("Can't allocated data\n");
return -1; 58 71 return -1;
} 59 72 }
60 73
printf ( "\n" ); 61 74 printf ( "\n" );
printf ( "HEATED_PLATE\n" ); 62 75 printf ( "HEATED_PLATE\n" );
printf ( " Parallel OpenMP version, using %d Threads\n",omp_get_max_threads() ); 63 76 printf ( " Parallel OpenMP version, using %d Threads\n",omp_get_max_threads() );
printf ( " A program to solve for the steady state temperature distribution\n" ); 64 77 printf ( " A program to solve for the steady state temperature distribution\n" );
printf ( " over a rectangular plate.\n" ); 65 78 printf ( " over a rectangular plate.\n" );
printf ( " Spatial grid of %d by %d points.\n\n", M, N ); 66 79 printf ( " Spatial grid of %d by %d points.\n\n", M, N );
67 80
68 81
/* Initialize grid and create input file */ 69 82 /* Initialize grid and create input file */
printf("Initializing grid\n"); 70 83 printf("Initializing grid\n");
71 84
inidat(N, M,u,unew); 72 85 inidat(N, M,u,unew);
73 86
prtdat(N, M,u, "initial.dat"); 74 87 prtdat(N, M,u, "initial.dat");
75 88
76
printf("Start computing\n"); 77 89 printf("Start computing\n");
78 90
int iter=0; 79 91 int iter=0;
80 92
/* 81 93 /*
* iterate until the new solution unew differs from the old solution u 82 94 * iterate until the new solution unew differs from the old solution u
* by no more than EPSILON. 83 95 * by no more than EPSILON.
* */ 84 96 * */
85 97
while(diff> EPSILON) { 86 98 while(diff> EPSILON) {
87 99
update(N, M, u, unew,&diff); 88 100 update(N, M, u, unew,&diff);
89 101
if(iter%ITER_PRINT==0) 90 102 if(iter%ITER_PRINT==0)
printf("Iteration %d, diff = %f\n ", iter,diff); 91 103 printf("Iteration %d, diff = %f\n ", iter,diff);
92 104
iter++; 93 105 iter++;
} 94 106 }
95 107
prtdat(N, M, u, "final.dat"); 96 108 prtdat(N, M, u, "final.dat");
97 109
free(u); 98 110 free(u);
free(unew); 99 111 free(unew);
} 100 112 }
101 113
102 114
103 115
/**************************************************************************** 104 116 /****************************************************************************
* subroutine update 105 117 * subroutine update
****************************************************************************/ 106 118 ****************************************************************************/
void update(int nx,int ny, float *u, float *unew, float * diff) 107 119 void update(int nx,int ny, float *u, float *unew, float * diff)
{ 108 120 {
int ix, iy; 109 121 int ix, iy;
*diff=0.0; 110 122 *diff=0.0;
111 123
#pragma omp parallel for shared(nx,ny,u,unew) private (ix,iy) 112 124 #pragma omp parallel for shared(nx,ny,u,unew) private (ix,iy)
for (ix = 1; ix < nx-1; ix++) { 113 125 for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) { 114 126 for (iy = 1; iy < ny-1; iy++) {
unew[ix*ny+iy] = 115 127 unew[ix*ny+iy] =
(u[(ix+1)*ny+iy] + u[(ix-1)*ny+iy] + 116 128 (u[(ix+1)*ny+iy] + u[(ix-1)*ny+iy] +
u[ix*ny+iy+1] + u[ix*ny+iy-1] )/4.0; 117 129 u[ix*ny+iy+1] + u[ix*ny+iy-1] )/4.0;
118 130
} 119 131 }
} 120 132 }
121 133
//compute reduction 122 134 //compute reduction
123 135
124 136
float mydiff; 125 137 float mydiff;
126 138
#pragma omp parallel shared(nx,ny,u,unew, diff) private (ix,iy,mydiff) 127 139 #pragma omp parallel shared(nx,ny,u,unew, diff) private (ix,iy,mydiff)
{ 128 140 {
mydiff=0.0; 129 141 mydiff=0.0;
#pragma omp for 130 142 #pragma omp for
for (ix = 1; ix < nx-1; ix++) { 131 143 for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) { 132 144 for (iy = 1; iy < ny-1; iy++) {
if (mydiff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] )) 133 145 if (mydiff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] ))
{ 134 146 {
mydiff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] ); 135 147 mydiff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] );
148 }
} 136 149 }
} 137 150 }
} 138
139 151
140 152
# pragma omp critical 141 153 # pragma omp critical
{ 142 154 {
if (*diff < mydiff ) 143 155 if (*diff < mydiff )
{ 144 156 {
*diff = mydiff; 145 157 *diff = mydiff;
} 146 158 }
} 147
148
149
#pragma omp for 150
for (ix = 1; ix < nx-1; ix++) { 151
for (iy = 1; iy < ny-1; iy++) { 152
u[ix*ny+iy] = unew[ix*ny+iy]; 153
} 154 159 }
} 155
156 160
157 161
158 162 #pragma omp for
159 163 for (ix = 1; ix < nx-1; ix++) {
} 160 164 for (iy = 1; iy < ny-1; iy++) {
165 u[ix*ny+iy] = unew[ix*ny+iy];
166 }
167 }
168 }
} 161 169 }
162 170
/***************************************************************************** 163 171 /*****************************************************************************
* Initialize Data 164 172 * Initialize Data
*****************************************************************************/ 165 173 *****************************************************************************/
void inidat(int nx, int ny, float *u, float *unew) 166 174 void inidat(int nx, int ny, float *u, float *unew)
{ 167 175 {
int ix, iy; 168 176 int ix, iy;
169 177
/* 170 178 /*
*Set boundary data and interrior values 171 179 *Set boundary data and interrior values
* */ 172 180 * */
for (ix = 0; ix < nx; ix++) 173 181 for (ix = 0; ix < nx; ix++)
for (iy = 0; iy < ny; iy++) { 174 182 for (iy = 0; iy < ny; iy++) {
175 183
if(ix==0) 176 184 if(ix==0)
{ 177 185 {
u[ix*ny+iy]=0.0; 178 186 u[ix*ny+iy]=0.0;
} 179 187 }
else 180 188 else
if(iy==0 && ix!=0) 181 189 if(iy==0 && ix!=0)
{ 182 190 {
u[ix*ny+iy]=100.0; 183 191 u[ix*ny+iy]=100.0;
}else 184 192 }else
185
if(ix==nx-1) 186
{ 187
u[ix*ny+iy]=100.0; 188
}else 189
190 193
if(iy==ny-1 && ix!=0) 191 194 if(ix==nx-1)
{ 192 195 {
u[ix*ny+iy]=100.0; 193 196 u[ix*ny+iy]=100.0;
}else 194 197 }else
195 198
u[ix*ny+iy]=( float ) ( 2 * nx + 2 * ny - 4 ); 196 199 if(iy==ny-1 && ix!=0)
200 {
201 u[ix*ny+iy]=100.0;
202 }else
203
204 u[ix*ny+iy]=0.0;
} 197 205 }
} 198 206 }
199 207
/************************************************************************** 200 208 /**************************************************************************
* Print Data to files 201 209 * Print Data to files
**************************************************************************/ 202 210 **************************************************************************/
void prtdat(int nx, int ny, float *u,const char *fnam) 203 211 void prtdat(int nx, int ny, float *u,const char *fnam)
{ 204 212 {
205 213
int ix, iy; 206 214 int ix, iy;
FILE *fp; 207 215 FILE *fp;
208 216
if(ITER_PRINT==0)return; 209 217 if(ITER_PRINT==0)return;
210 218
fp = fopen(fnam, "w"); 211 219 fp = fopen(fnam, "w");
212 220
for (ix = 0 ; ix < nx; ix++) { 213 221 for (ix = 0 ; ix < nx; ix++) {
for (iy =0; iy < ny; iy++) { 214 222 for (iy =0; iy < ny; iy++) {
lab3/ser_heat2D.c View file @ d16e53a
/**************************************************************************** 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 +------------------+
u = 100 20 20 u = 100
21 21
Interrior point : 22 22 Interrior point :
u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] ) 23 23 u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] )
24 24
****************************************************************************/ 25 25 ****************************************************************************/
#include <stdio.h> 26 26 #include <stdio.h>
#include <stdlib.h> 27 27 #include <stdlib.h>
#include <math.h> 28 28 #include <math.h>
29 29
#define N 200 30 30 #define NN 50
#define M 200 31 31 #define MM 50
32 32
#define ITER_PRINT 100 33 33 #define ITER_PRINT 100
#define PRINT_DATA 1 34 34 #define PRINT_DATA 1
35 35
#define EPSILON 1e-1 36 36 #define _EPSILON 0.01
37 37
38 38
void update(int nx,int ny, float *u, float *unew, float * diff); 39 39 void update(int nx,int ny, float *u, float *unew, float * diff);
void inidat(int nx, int ny, float *u, float *unew); 40 40 void inidat(int nx, int ny, float *u, float *unew);
void prtdat(int nx, int ny, float *u,const char *fnam); 41 41 void prtdat(int nx, int ny, float *u,const char *fnam);
42 42
43 43
44 44
45 45
int main(int argc, char *argv[]) 46 46 int main(int argc, char *argv[])
{ 47 47 {
48 int N=NN,M=MM;
48 49
50 float EPSILON=_EPSILON;
51
52 if(argc !=3)
53 {
54 fprintf(stderr,"usage %s N EPSILON\n ", argv[0]);
55 fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n");
56 fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n");
57 return -1;
58 }
59
60 N = M = atoi(argv[1]);
61 EPSILON = atof(argv[2]);
62
float diff=1.0; 49 63 float diff=1.0;
50 64
float *u = (float *)malloc(N * M * sizeof(float)); 51 65 float *u = (float *)malloc(N * M * sizeof(float));
float *unew = (float *)malloc(N * M * sizeof(float)); 52 66 float *unew = (float *)malloc(N * M * sizeof(float));
53 67
if(u==0 || unew ==0) 54 68 if(u==0 || unew ==0)
{ 55 69 {
perror("Can't allocated data\n"); 56 70 perror("Can't allocated data\n");
return -1; 57 71 return -1;
} 58 72 }
59 73
printf ( "\n" ); 60 74 printf ( "\n" );
printf ( "HEATED_PLATE\n" ); 61 75 printf ( "HEATED_PLATE\n" );
printf ( " Serial version\n" ); 62 76 printf ( " Serial version\n" );
printf ( " A program to solve for the steady state temperature distribution\n" ); 63 77 printf ( " A program to solve for the steady state temperature distribution\n" );
printf ( " over a rectangular plate.\n" ); 64 78 printf ( " over a rectangular plate.\n" );
printf ( "\n" ); 65 79 printf ( "\n" );
printf ( " Spatial grid of %d by %d points.\n\n", M, N ); 66 80 printf ( " Spatial grid of %d by %d points.\n", M, N );
81 printf ( " The iteration will end until tolerance <= %f\n\n",EPSILON);
67 82
68
/* Initialize grid and create input file */ 69 83 /* Initialize grid and create input file */
printf("Initializing grid\n"); 70 84 printf("Initializing grid\n");
71 85
inidat(N, M,u,unew); 72 86 inidat(N, M,u,unew);
73 87
prtdat(N, M,u, "initial.dat"); 74 88 prtdat(N, M,u, "initial.dat");
75
76 89
printf("Start computing\n"); 77 90 printf("Start computing\n\n");
78 91
int iter=0; 79 92 int iter=0;
80 93
/* 81 94 /*
* iterate until the new solution unew differs from the old solution u 82 95 * iterate until the new solution unew differs from the old solution u
* by no more than EPSILON. 83 96 * by no more than EPSILON.
* */ 84 97 * */
85 98
while(diff> EPSILON) { 86 99 while(diff> EPSILON) {
87 100
update(N, M, u, unew,&diff); 88 101 update(N, M, u, unew,&diff);
89 102
if(iter%ITER_PRINT==0) 90 103 if(iter%ITER_PRINT==0)
printf("Iteration %d, diff = %f\n ", iter,diff); 91 104
105 printf("Iteration %d, diff = %f\n ", iter,diff);
92 106
iter++; 93 107 iter++;
} 94 108 }
95 109
prtdat(N, M, u, "final.dat"); 96 110 prtdat(N, M, u, "final.dat");
97 111
free(u); 98 112 free(u);
free(unew); 99 113 free(unew);
} 100 114 }
101 115
102 116
103 117
/**************************************************************************** 104 118 /****************************************************************************
* subroutine update 105 119 * subroutine update
****************************************************************************/ 106 120 ****************************************************************************/
void update(int nx,int ny, float *u, float *unew, float * diff) 107 121 void update(int nx,int ny, float *u, float *unew, float * diff)
{ 108 122 {
int ix, iy; 109 123 int ix, iy;
*diff=0.0; 110 124 *diff=0.0;
111 125
for (ix = 1; ix < nx-1; ix++) { 112 126 for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) { 113 127 for (iy = 1; iy < ny-1; iy++) {
unew[ix*ny+iy] = 114 128 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
(u[(ix+1)*ny+iy] + u[(ix-1)*ny+iy] + 115 129 ;
u[ix*ny+iy+1] + u[ix*ny+iy-1] )/4.0; 116
117
if (*diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] )) 118 130 if (*diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] ))
{ 119 131 {
*diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] ); 120 132 *diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] );
} 121 133 }
} 122 134 }
123 135
} 124 136 }
125 137
126 138
for (ix = 1; ix < nx-1; ix++) { 127 139 for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) { 128 140 for (iy = 1; iy < ny-1; iy++) {
u[ix*ny+iy] = unew[ix*ny+iy]; 129 141 u[ix*ny+iy] = unew[ix*ny+iy];
} 130 142 }
} 131 143 }
132 144
} 133 145 }
134 146
/***************************************************************************** 135 147 /*****************************************************************************
* Initialize Data 136 148 * Initialize Data
*****************************************************************************/ 137 149 *****************************************************************************/
void inidat(int nx, int ny, float *u, float *unew) 138 150 void inidat(int nx, int ny, float *u, float *unew)
{ 139 151 {
int ix, iy; 140 152 int ix, iy;
141 153
/* 142 154 /*
*Set boundary data and interrior values 143 155 *Set boundary data and interrior values
* */ 144 156 * */
for (ix = 0; ix < nx; ix++) 145 157 for (ix = 0; ix < nx; ix++)
for (iy = 0; iy < ny; iy++) { 146 158 for (iy = 0; iy < ny; iy++) {
147 159
if(ix==0) 148 160 if(ix==0)
{ 149 161 {
u[ix*ny+iy]=0.0; 150 162 u[ix*ny+iy]=0.0;
} 151 163 }
else 152 164 else
if(iy==0 && ix!=0) 153 165 if(iy==0 && ix!=0)
{ 154 166 {
u[ix*ny+iy]=100.0; 155 167 u[ix*ny+iy]=100.0;
}else 156 168 }else
157 169
if(ix==nx-1) 158 170 if(ix==nx-1)
{ 159 171 {
u[ix*ny+iy]=100.0; 160 172 u[ix*ny+iy]=100.0;
}else 161 173 }else
162 174
if(iy==ny-1 && ix!=0) 163 175 if(iy==ny-1 && ix!=0)
{ 164 176 {
u[ix*ny+iy]=100.0; 165 177 u[ix*ny+iy]=100.0;
}else 166 178 }else
167 179
u[ix*ny+iy]=( float ) ( 2 * nx + 2 * ny - 4 ); 168 180 u[ix*ny+iy]=0.0;