source: CIVL/examples/translation/mpi/diffusion2d.c@ 1778e7a

1.23 2.0 acw/focus-triggers main test-branch
Last change on this file since 1778e7a was 696be79, checked in by Ziqing Luo <ziqing@…>, 12 years ago

changed diffusion2d example filename," xProcs " to "NPROCSX" ditto for "y". Get rid of all #ifdefs in function setBoundaries().

Corresponding changes for Makefile and tests

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@1745 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 10.0 KB
Line 
1/* diffusion2d_cb.c: parallel 2d-diffusion equation solver with constant boundaries
2 * slicing matrix as a checker board.
3 * To execute: mpicc diffusion2d_cb.c ; mpiexec -n 4 ./a.out 2 2
4 * To verify: civl verify diffusion2d_cb.c
5 */
6#include<stdio.h>
7#include<stdlib.h>
8#include<assert.h>
9#include<string.h>
10#include<mpi.h>
11
12/* Message tags */
13#define FROMLEFT 0
14#define FROMRIGHT 1
15#define FROMTOP 2
16#define FROMBOTTOM 3
17#define DATAPASS 4
18#define comm MPI_COMM_WORLD
19
20#ifdef _CIVL
21
22$input int NXB = 5; // nx upper bound
23$input int nx; // global number of columns in matrix
24$assume 1 <= nx && nx <= NXB;
25$input int NYB = 5; // ny upper bound
26$input int ny; // global number of rows of matrix
27$assume 1 <= ny && ny <= NYB;
28$input double u_init[ny+2][nx+2]; // initial value of temperatures, includes 4
29 // strings of constant boundries
30$input double k; // constant coefficient
31$assume k > 0.0 && k < 0.5;
32$input int NSTEPSB = 5; // boundary of number of steps
33$input int nsteps; // number of steps
34$assume 1<=nsteps && nsteps<=NSTEPSB;
35$input int wstep = 1; // write frame every this many time steps
36double oracle[nsteps][ny+2][nx+2]; // solution computed sequentially, done by proc 0 only
37$input int XPROCSB; // Bound number of components of columns
38$input int NPROCSX; // Number of components of columns
39$assume NPROCSX > 1 && NPROCSX <= XPROCSB;
40$input int YPROCSB; // Bound number of components of rows
41$input int NPROCSY; // Number of components of rows
42$assume NPROCSY > 1 && NPROCSY <= YPROCSB;
43$input int _NPROCS;
44#else
45
46int nx, ny, nsteps, wstep;
47int NPROCSX, NPROCSY;
48double constTemp, initTemp; // values of constant boundaries and
49 // initial tempretures
50double k;
51
52#endif
53
54/* Global variables */
55double ** u_curr;
56double ** u_next;
57int nprocs, rank, left, right, top, bottom;
58int nxl, nyl, firstCol, firstRow;
59
60/* Compute the global column index of cells owned by the process */
61int firstColForProc(int rank) {
62 return (rank - (rank / NPROCSX)*NPROCSX)*nx/NPROCSX;
63}
64
65/* Compute the global row index of cells owned by the process */
66int firstRowForProc(int rank) {
67 return ((rank / NPROCSX)*ny)/NPROCSY;
68}
69
70/* Computes the number of columns owned by the process */
71int countColForProc(int rank) {
72 int a = firstColForProc(rank);
73 int b;
74
75 if((rank / NPROCSX) == ((rank+1) / NPROCSX))
76 b = firstColForProc(rank+1);
77 else
78 b = nx;
79 return b - a;
80}
81
82/* Computes the number of rows owned by the process */
83int countRowForProc(int rank) {
84 int a = firstRowForProc(rank);
85 int b = firstRowForProc(rank+NPROCSX);
86
87 return b - a;
88}
89
90/* Get the owner process of the given cell */
91int OWNER(int col, int row) {
92 int procRow = ((NPROCSY * (row+1))-1) / ny;
93 int procCol = ((col + 1) * NPROCSX - 1) / nx;
94
95 return procRow * NPROCSX + procCol;
96}
97
98
99#ifdef _CIVL
100void setConstBoundaries() {
101 int i;
102
103 // sets vertical constant boundaries
104 if(left == MPI_PROC_NULL)
105 for(i=0; i<nyl+2; i++) {
106 u_curr[i][0] = u_init[i + firstRow][0];
107 u_next[i][0] = u_curr[i][0];
108 }
109 if(right == MPI_PROC_NULL)
110 for(i=0; i<nyl+2; i++) {
111 u_curr[i][nxl+1] = u_init[i + firstRow][nx+1];
112 u_next[i][nxl+1] = u_curr[i][nxl+1];
113 }
114 // sets horizontal constant boundaries
115 if(top == MPI_PROC_NULL)
116 for(i=0; i<nxl+2; i++) {
117 u_curr[0][i] = u_init[0][i + firstCol];
118 u_next[0][i] = u_curr[0][i];
119 }
120 if(bottom == MPI_PROC_NULL)
121 for(i=0; i<nxl+2; i++) {
122 u_curr[nyl+1][i] = u_init[ny+1][i + firstCol];
123 u_next[nyl+1][i] = u_curr[nyl+1][i];
124 }
125}
126#endif
127
128/* Initialize all global variables */
129void initialization(int argc, char * argv[]) {
130 int i,j;
131
132#ifndef _CIVL
133
134 nsteps = 300;
135 wstep = 5;
136 nx = 15;
137 ny = 15;
138 if(argc < 3) {
139 printf("Program needs 2 arguments to specify the number of sliced components in x axis and y axis:\n"
140 "It should go with the format: mpiexec -n [nprocs] "
141 "[filename] [#components in x axis] [#components in y axis]\n");
142 assert(0);
143 }
144 NPROCSX = atoi(argv[1]);
145 NPROCSY = atoi(argv[2]);
146 assert(0 < NPROCSX * NPROCSY <=nprocs);
147 nprocs = (nprocs > NPROCSX * NPROCSY)?NPROCSX*NPROCSY:nprocs;
148 constTemp = 0.0;
149 initTemp = 100.0;
150 k = 0.13;
151 printf("Diffusion2d with k=%f, nx=%d, ny=%d, nsteps=%d, wstep=%d\n",
152 k, nx, ny, nsteps, wstep);
153
154#endif
155
156 nxl = countColForProc(rank);
157 nyl = countRowForProc(rank);
158 u_curr = (double **)malloc((nyl + 2) * sizeof(double *));
159 assert(u_curr);
160 u_next = (double **)malloc((nyl + 2) * sizeof(double *));
161 assert(u_next);
162 for(i=0; i < nyl + 2; i++){
163 u_curr[i] = (double *)malloc((nxl + 2) * sizeof(double));
164 assert(u_curr[i]);
165 u_next[i] = (double *)malloc((nxl + 2) * sizeof(double));
166 assert(u_next[i]);
167 }
168 firstCol = firstColForProc(rank);
169 firstRow = firstRowForProc(rank);
170 // computes neighbors
171 if(firstCol != 0)
172 left = OWNER(firstCol - 1, firstRow);
173 else
174 left = MPI_PROC_NULL;
175 if(firstRow != 0)
176 top = OWNER(firstCol, firstRow - 1);
177 else
178 top = MPI_PROC_NULL;
179 if(firstCol + nxl < nx)
180 right = OWNER(firstCol + nxl, firstRow);
181 else
182 right = MPI_PROC_NULL;
183 if(firstRow + nyl < ny)
184 bottom = OWNER(firstCol, firstRow + nyl);
185 else
186 bottom = MPI_PROC_NULL;
187#ifdef _CIVL
188 setConstBoundaries();
189 // In CIVL mode process with rank 0 will be responsible for computing the diffusion2d equation
190 // sequentially such that the results can be used to compare with the ones of parallel run.
191 if(rank == 0) {
192 for(i = 0; i < ny + 2; i++)
193 for(j = 0; j < nx + 2; j++)
194 oracle[0][i][j] = u_init[i][j];
195 for(int t=1; t < nsteps; t++)
196 for(i = 0; i < ny + 2; i++)
197 for(j = 0; j < nx + 2; j++)
198 if(i==0 || j==0 || i == ny + 1 || j == nx + 1)
199 oracle[t][i][j] = oracle[t-1][i][j];
200 else
201 oracle[t][i][j] = oracle[t-1][i][j] +
202 k*(oracle[t-1][i+1][j] + oracle[t-1][i-1][j] +
203 oracle[t-1][i][j+1] + oracle[t-1][i][j-1] - 4*oracle[t-1][i][j]);
204 }
205
206#endif
207}
208
209void update() {
210 int i, j;
211 double **tmp;
212
213 for(i = 1; i < nyl + 1; i++)
214 for(j = 1; j < nxl + 1; j++) {
215 u_next[i][j] = u_curr[i][j] +
216 k*(u_curr[i+1][j] + u_curr[i-1][j] +
217 u_curr[i][j+1] + u_curr[i][j-1] - 4*u_curr[i][j]);
218 }
219 //swap two pointers
220 tmp = u_curr;
221 u_curr = u_next;
222 u_next = tmp;
223}
224
225void exchange() {
226 double sendbuf[nyl];
227 double recvbuf[nyl];
228
229 // sends top string, receives bottom string
230 MPI_Sendrecv(&u_curr[1][1], nxl, MPI_DOUBLE, top, FROMBOTTOM, &u_curr[nyl+1][1], nxl,
231 MPI_DOUBLE, bottom, FROMBOTTOM, comm, MPI_STATUS_IGNORE);
232 // sends bottom string, receives top string
233 MPI_Sendrecv(&u_curr[nyl][1], nxl, MPI_DOUBLE, bottom, FROMTOP, &u_curr[0][1], nxl,
234 MPI_DOUBLE, top, FROMTOP, comm, MPI_STATUS_IGNORE);
235 // sends left most string, receives right most string
236 for(int i = 0; i < nyl; i++) sendbuf[i] = u_curr[i+1][1];
237 MPI_Sendrecv(sendbuf, nyl, MPI_DOUBLE, left, FROMRIGHT, recvbuf, nyl,
238 MPI_DOUBLE, right, FROMRIGHT, comm, MPI_STATUS_IGNORE);
239 if(right != MPI_PROC_NULL)
240 for(int i = 0; i < nyl; i++) u_curr[i+1][nxl+1] = recvbuf[i];
241 // sends right most string, receives left most string
242 for(int i = 0; i < nyl; i++) sendbuf[i] = u_curr[i+1][nxl];
243 MPI_Sendrecv(sendbuf, nyl, MPI_DOUBLE, right, FROMLEFT, recvbuf, nyl,
244 MPI_DOUBLE, left, FROMLEFT, comm, MPI_STATUS_IGNORE);
245 if(left != MPI_PROC_NULL)
246 for(int i = 0; i < nyl; i++) u_curr[i+1][0] = recvbuf[i];
247}
248
249void printData(int time, int firstCol, int nxl, int firstRow, int nyl, double * buf) {
250
251 for(int i=0; i<nyl; i++) {
252 for(int j=0; j<nxl; j++) {
253 printf("%8.2f ", *(buf + i*nxl + j));
254#ifdef _CIVL
255 $assert(*(buf + i*nxl + j) == oracle[time][firstRow + i + 1][firstCol + j + 1]) : \
256 "Error: disagreement at time %d position [%d][%d]: saw %lf, expected %lf", \
257 time, firstRow + i, firstCol + j,
258 *(buf + i*nxl + j), oracle[time][firstRow + i + 1][firstCol + j + 1];
259#endif
260 }
261 printf("\n");
262 }
263}
264
265void write_frame(int time) {
266 double * buf; // buffer of data to print
267 int i, j;
268
269 buf = (double *)malloc(nxl * nyl * sizeof(double));
270 // writes data into buffer array
271 for(i = 0; i < nyl; i++)
272 for(j = 0; j < nxl; j++) {
273 buf[i*nxl + j] = u_curr[i+1][j+1];
274 }
275 if(rank == 0) {
276 printf("\n-------------------- time step:%d --------------------\n", time);
277 printData(time, firstCol, nxl, firstRow, nyl, buf);
278 free(buf);
279 for(i=1; i<nprocs; i++){
280 double * recvbuf;
281 int senderx, sendery;
282 int senderNyl, senderNxl;
283
284 senderNxl = countColForProc(i);
285 senderNyl = countRowForProc(i);
286 if(senderNxl != 0 && senderNyl != 0) {
287 recvbuf = (double *)malloc(senderNxl * senderNyl * sizeof(double));
288 senderx = firstColForProc(i);
289 sendery = firstRowForProc(i);
290 MPI_Recv(recvbuf, senderNyl*senderNxl, MPI_DOUBLE, i,
291 DATAPASS, comm, MPI_STATUS_IGNORE);
292 printData(time, senderx, senderNxl, sendery, senderNyl, recvbuf);
293 free(recvbuf);
294 }
295 }
296 } else {
297 MPI_Send(buf, nyl*nxl, MPI_DOUBLE, 0, DATAPASS, comm);
298 free(buf);
299 }
300}
301
302int main(int argc, char * argv[]) {
303 int i,j;
304
305#ifdef _CIVL
306
307 // elaborating nx, ny, NPROCSX and NPROCSY...
308 elaborate(nx);
309 elaborate(ny);
310 elaborate(NPROCSX);
311 elaborate(NPROCSY);
312
313#endif
314 MPI_Init(&argc, &argv);
315 MPI_Comm_rank(comm, &rank);
316 MPI_Comm_size(comm, &nprocs);
317 initialization(argc, argv);
318#ifdef _CIVL
319
320 for(i=1; i<nyl+1; i++)
321 memcpy(&u_curr[i][1], &u_init[firstRow + i][firstCol + 1], nxl * sizeof(double));
322
323#else
324
325 for(int i=0; i < nyl+2; i++)
326 for(int j=0; j < nxl+2; j++)
327 if(i == 0 || j == 0 || i == nyl+1 || j == nxl+1)
328 u_next[i][j] = u_curr[i][j] = constTemp;
329 else
330 u_curr[i][j] = initTemp;
331
332#endif
333 for(i=0; i<nsteps; i++) {
334 if(nxl != 0 && nyl != 0) {
335 if(i%wstep == 0)
336 write_frame(i);
337 exchange();
338 update();
339 }
340 }
341 for(i=0; i<nyl+2; i++) {
342 free(u_curr[i]);
343 free(u_next[i]);
344 }
345 free(u_curr);
346 free(u_next);
347 return 0;
348}
Note: See TracBrowser for help on using the repository browser.