Category Archives: MPI

MPI/C – 2 Dimensional Heat Equation [Square Grid]

MPI

For my understanding of what MPI is &/or does, please refer to this post.

Program Listing

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
/* heat_equation_2d.c
 * PARALLEL [MPI] C PROGRAM TO SOLVE THE 2 DIMENSIONAL HEAT EQUATION.
 * THE INITIAL TEMPERATURE IS COMPUTED TO BE HIGH IN THE MIDDLE OF THE
 * DOMAIN UNDER CONSIDERATION AND ZERO AT THE BOUNDARIES. THE BOUNDARIES
 * ARE HELD AT ZERO TEMPERATURE THROUGH OUT THE SIMULATION.
 *
 * THE GIVEN GRID IS DECOMPOSED BY THE MASTER AND THEN DISTRIBUTED BY
 * COLUMNS TO WORKERS. A GRID POINT'S CURRENT TEMPERATURE DEPENDS UPON 
 * ITS PREVIOUS TIME STEP VALUE AS WELL AS THE VALUES OF THE NEIGHBORING
 * GRID POINTS. AS SUCH, AT EACH TIME STEP, WORKERS EXCHANGE BORDER/BOUNDARY
 * DATA WITH NEIGHBORS. ONCE ALL THE TIME STEPS ARE COMPLETED, WORKERS
 * SEND THEIR RESULTS TO MASTER.
 *
 * THE PROGRAM GENERATES A DATA FILE THAT CAN BE VISUALIZED WITH THIRD 
 * PARTY PROGRAMS [GNUPLOT, MATLAB, ETC.]
 *
 * TESTED SUCCESSFULLY WITH MPICH2 (1.3.1) COMPILED AGAINST GCC (4.1.2) 
 * IN A LINUX BOX WITH QUAD CORE INTEL XEON PROCESSOR (1.86 GHz) & 4GB OF RAM.
 *
 * ORIGINAL AUTHOR : BLAISE BARNEY; 2009/12/11
 *                   LAWRENCE LIVERMORE NATIONAL LABORATORY
 * FIRST COPIED    : GOWTHAM; Fri, 26 Nov 2010 21:30:30 -0500
 * LAST MODIFIED   : GOWTHAM; Fri, 10 Dec 2010 01:34:52 -0500
 *
 * URL:
 * http://sgowtham.com/journal/2010/12/11/mpi-c-2-dimensional-heat-equation/
 *
 * COMPILATION:
 * mpicc -g -Wall -lm heat_equation_2d.c -o heat_equation_2d.x 
 *
 * EXECUTION:
 * mpirun -np NPROC -machinefile MACHINEFILE ./heat_equation_2d.x
 *
 * NPROC       : NUMBER OF PROCESSORS ALLOCATED TO RUNNING THIS PROGRAM;
 *               NXY [GRID SIZE] MUST BE DIVISIBLE BY [NPROC - 1]
 * MACHINEFILE : FILE LISTING THE HOSTNAMES OF PROCESSORS ALLOCATED TO
 *               RUNNING THIS PROGRAM
 *
 * OUTPUT      : DATA FILE 'heat_equation_2d.dat' WITH THE FOLLOWING
 *               FORMAT
 *
 *               TSTEPS=#TSTEPS
 *               U(x, y) AT t = 0         [NXY x NXY]
 *               U(x, y) AT t = 1         [NXY x NXY]
 *               ...
 *               U(x, y) AT t = TSTEP     [NXY x NXY]
 *
 *               TOTAL NUMBER OF LINES IN 'heat_equation_2d.dat'
 *               (NXY * (TSTEPS + 1)) + 1
 *
*/
 
/* STANDARD HEADERS AND DEFINITIONS 
 * REFERENCE: http://en.wikipedia.org/wiki/C_standard_library
*/
#include <stdio.h>  /* Core input/output operations                         */
#include <stdlib.h> /* Conversions, random numbers, memory allocation, etc. */
#include <math.h>   /* Common mathematical functions                        */
#include <time.h>   /* Converting between various date/time formats         */
#include <mpi.h>    /* MPI functionality                                    */
 
#define NXY          20          /* Number of (x or y) grids in U(x,y) */
                                 /* NXY must be divisible by NPROC - 1 */
#define TSTEPS      500          /* Number of time steps               */
#define MASTER        0          /* Process ID for MASTER              */
#define BEGIN         1          /* Message tag                        */
#define DONE          2          /* Message tag                        */
#define LNEIGHBOR     3          /* Left neighbor tag                  */
#define RNEIGHBOR     4          /* Right neigbor tag                  */
#define NONE          5          /* Indicates no neighbor              */
 
struct ThermalDiffusivity {      /* Parameters used in difference      */
  double cx;                     /* equations - refer to class         */
  double cy;                     /* notes. Physically, these refer to  */
} parameters = {0.25, 0.25};     /* thermal diffusivity                */
 
 
/* MAIN PROGRAM BEGINS */
int main (int argc, char *argv[]) {
 
  /* VARIABLE DECLARATION & INITIALIZATION */
  double U[NXY][NXY],      /* 2D array to hold current values        */
         Unew[NXY][NXY],   /* 2D array to hold new values            */
         start_time,       /* Wall clock - start time                */
         end_time;         /* Wall clock - end time                  */
 
  int    proc_id,          /* Process ID                             */
         n_procs,          /* Number of processors                   */
         n_workers,        /* Number of WORKERS [nprocs - 1]         */
         cols_per_worker,  /* Number of columns of U(x,y) per WORKER */ 
         offset,           /* Offset                                 */
         source,           /* Source for a message                   */
         destination,      /* Destination for a message              */
         lneighbor,        /* Process ID of left neighbor            */
         rneighbor,        /* Process ID of right neighbor           */
         start,            /* Starting # for columns in WORKERS      */
         end,              /* Ending # for columns in WORKERS        */
         i,                /* Dummy/Running index                    */
         j,                /* Dummy/Running index                    */
         h;                /* Dummy/Running index for TSTEPS         */
 
  FILE   *output;          /* File pointer                           */
 
  MPI_Status status;       /* MPI structure containing return codes
                              for message passing operations         */
 
  MPI_Datatype cvector;    /* A custom MPI data type defined to 
                              contain column vectors                 */
 
 
  /* INITIALIZE MPI */
  MPI_Init(&argc, &argv);
 
  /* GET THE PROCESSOR ID AND NUMBER OF PROCESSORS */
  MPI_Comm_rank(MPI_COMM_WORLD, &proc_id);
  MPI_Comm_size(MPI_COMM_WORLD, &n_procs);
 
  /* NUMBER OF WORKERS [NPROC - 1] */
  n_workers = n_procs - 1;
 
  /* MPI DATA TYPE [CUSTOM] - TO SEND BLOCKS OF U(x,y) TO WORKERS AS COLUMNS
   * THANKS TO DR. STEVEN SEIDEL [http://www.cs.mtu.edu/~steve/] 
   * FOR INFORMING/TEACHING ME THE POSSIBILITY OF DEFINING CUSTOM DATA TYPES 
   * IN MPI, DURING A PREVIOUS DISCUSSION ABOUT MATRIX MULTIPILICATION
  */
  MPI_Type_vector(NXY, 1, NXY, MPI_DOUBLE, &cvector);
  MPI_Type_commit (&cvector);
 
  /* IF MASTER, THEN DO THE FOLLOWING:
   * PRINT GENERAL INFORMATION ABOUT THE PROGRAM
   * CHECK IF NPROCS MEETS THE PROGRAM REQUIREMENTS
   * POPULATE THE U(x,y) TO SIMULATE INITIAL TEMPERATURE DISTRIBUTION
   * SEND RELEVANT INFORMATION ABOUT NEIGHBORS, OF U(x,y) TO WORKERS
   * RECEIVE RELAVANT INFORMATION FROM WORKERS, AT EVERY TIME STEP
   * DISPLAY U(x,y) AT INITIAL AND FINAL STEPS
   * PRINT U(x,y) AT EVERY TIME STEP SO THAT IT CAN BE PLOTTED/ANIMATED
  */
  if (proc_id == MASTER) {
 
    printf("\n  2D Heat Equation [MPI/C]\n\n");
    printf("    Number of x-grids        : %d\n", NXY);
    printf("    Number of y-grids        : %d\n", NXY);
    printf("    Number of processors     : %d\n", n_procs);
    printf("    Number of workers        : %d\n", n_workers);
    printf("    Number of time steps     : %d\n", TSTEPS);
    printf("    Thermal diffusivity (Cx) : %4.3lf\n", parameters.cx);
    printf("    Thermal diffusivity (Cy) : %4.3lf\n\n", parameters.cy);
 
    /* CHECK IF NXY IS DIVISIBLE BY [NPROC + 1]
     * IF NOT, PRINT AN ERROR MESSAGE AND QUIT THE PROGRAM 
    */
    if (NXY%(n_procs - 1) != 0) {
      printf("    ERROR:\n");
      printf("    NXY not divisible by (NPROC + 1)\n\n");
 
      exit(1);
    }
 
    /* OPEN THE FILE FOR WRITING [OVER-WRITE, IF EXISTS]
     * PRINT ERROR MESSAGE, IF ANY & TERMINATE THE PROGRAM
    */
    if(!(output = fopen("heat_equation_2d.dat", "w"))) {
      fprintf(stderr, "\n\n    ERROR:\n");
      fprintf(stderr, "    Unable to open the output file\n\n");
 
      exit(1);
    }
 
    /* START THE CLOCK */
    start_time = MPI_Wtime();
 
    /* INITIALIZE THE GRID */
    printf("    Temperature distribution (at time t=0)\n\n");
 
    /* HEADER LINE INDICATING TOTAL NUMBER OF TIME STEPS
     * NEEDED BY MATLAB (&/OR OTHER) FOR PLOTTING
    */ 
    fprintf(output, "TSTEPS=%d\n", TSTEPS + 1);
 
    /* POPULATE THE ARRAY, U(x,y) 
     * BORDER ELEMENTS ARE ALL ZERO & HIGH VALUE IN THE MIDDLE
     * PRINT THE ARRAY - TO THE SCREEN AND TO THE FILE
    */
    for (i = 0; i < NXY; i++) {
      for (j = 0; j < NXY; j++) {
        U[i][j] = 100 * i * (NXY - i - 1) * j * (NXY - j - 1);
        printf("    %08.0lf", U[i][j]);
        fprintf(output, "%012.4lf  ", U[i][j]);
      }
      printf("\n");
      fprintf(output, "\n");
    }
 
    /* NUMBER OF COLUMNS PER WORKER 
     * IT'S IMPORTANT TO ALLOCATE EQUAL NUMBER OF COLUMNS TO
     * EACH WORKER - ELSE, SOME MIGHT FINISH FASTER WITHOUT
     * ADDITIONAL WORK
     */
    cols_per_worker = NXY/n_workers;
 
    /* INITIALIZE offset */
    offset = 0;
 
    for (i = 1; i <= n_workers; i++) {
 
      /* INFORM EACH WORKER OF ITS NEIGHBORS */
      if (i == 1) {
        lneighbor = NONE;
      } else {
        lneighbor = i - 1;
      }
 
      if (i == n_workers) {
        rneighbor = NONE;
      } else {
        rneighbor = i + 1;
      }
 
      /* DESTINATION PROCESSOR */
      destination = i;
 
      /* SEND INFORMATION TO EACH WORKER
       * OFFSET
       * NUMBER OF COLUMNS EACH WORKER HAS TO WORK WITH 
       * LEFT & RIGHT NEIGHBORS FOR EACH WORKER
       * RELEVANT PORTION OF U(x,y) TO EACH WORKER
      */
      MPI_Send(&offset,    1, MPI_INT, destination, BEGIN, MPI_COMM_WORLD);
      MPI_Send(&cols_per_worker, 1, MPI_INT, destination, BEGIN, MPI_COMM_WORLD);
      MPI_Send(&lneighbor, 1, MPI_INT, destination, BEGIN, MPI_COMM_WORLD);
      MPI_Send(&rneighbor, 1, MPI_INT, destination, BEGIN, MPI_COMM_WORLD);
 
      for(j=0; j < cols_per_worker; j++) {
        MPI_Send(&U[0][offset+j], 1, cvector, destination, BEGIN, MPI_COMM_WORLD);
      }
 
      /* INCREMENT THE offset */
      offset = offset + cols_per_worker;
    }
 
    /* RECEIVE RELEVANT INFORMATION FROM WORKERS  FOR EVERY TIME STEP */
    for(h = 1; h <= TSTEPS; h++) {
 
      for(i = 1; i <= n_workers; i++) {
 
        /* SOURCE PROCESSOR */
        source = i;
 
        MPI_Recv(&offset, 1, MPI_INT, source, DONE, MPI_COMM_WORLD, &status); 
        MPI_Recv(&cols_per_worker, 1, MPI_INT, source, DONE, MPI_COMM_WORLD, &status); 
 
        for(j=0; j < cols_per_worker; j++) {
          MPI_Recv(&U[0][offset+j], 1, cvector, source, h*DONE, MPI_COMM_WORLD, &status);
        }
 
      }
 
      /* PRINT U(x,y), THE TEMPERATURE DISTRIBUTION AFTER EVERY TIME STEP TO THE FILE */
      for (i = 0; i < NXY; i++) {
        for (j = 0; j < NXY; j++) {
          fprintf(output, "%012.4lf  ", U[i][j]);
        }
        fprintf(output, "\n");
      }
 
    }
 
    /* CLOSE THE FILE */
    fclose(output);
 
    /* PRINT U(x,y), THE TEMPERATURE DISTRIBUTION AFTER 'TSTEPS' STEPS TO THE SCREEN */
    printf("\n\n    Temperature distribution (after %d time steps)\n\n", h-1);
    for (i = 0; i < NXY; i++) {
      for (j = 0; j < NXY; j++) {
        printf("    %08.0lf", U[i][j]);
      }
      printf("\n");
    }
 
    printf("\n\n");
 
    /* STOP THE CLOCK */
    end_time = MPI_Wtime();
 
    /* PRINT TIMING INFORMATION */
    printf("    Time elapsed (sec)   : %6.4f\n\n", end_time - start_time);
 
  } /* MASTER LOOP ENDS */
 
 
  /* IF WORKER, THEN DO THE FOLLOWING:
   * INITIALIZE U(x,y) & Unew(x,y) TO ZERO
   * RECEIVE RELEVANT INFORMATION FROM MASTER
   * CALCULATE U(x, y, t+1) BY COMMUNICATING APPROPRIATELY WITH NEIGHBORING WORKERS
   * SEND RELEVANT INFORMATION TO MASTER 
  */
  if (proc_id != MASTER) {
 
    /* INITIALIZE ALL ELEMENTS OF U(x,y) & Unew(x,y) TO ZERO */
    for(i = 0; i < NXY; i++) {
      for(j = 0; j < NXY; j++) {
        U[i][j] = 0.0;
        Unew[i][j] = 0.0;
      }
    }
 
    /* RECEIVE THE COLUMNS, NEIGHBORS AND GRID INFORMATION FROM MASTER */
    MPI_Recv(&offset,    1, MPI_INT, MASTER, BEGIN, MPI_COMM_WORLD, &status);
    MPI_Recv(&cols_per_worker, 1, MPI_INT, MASTER, BEGIN, MPI_COMM_WORLD, &status);
    MPI_Recv(&lneighbor, 1, MPI_INT, MASTER, BEGIN, MPI_COMM_WORLD, &status);
    MPI_Recv(&rneighbor, 1, MPI_INT, MASTER, BEGIN, MPI_COMM_WORLD, &status);
 
    for(j=0; j < cols_per_worker; j++) {
      MPI_Recv(&U[0][offset+j], 1, cvector, MASTER, BEGIN, MPI_COMM_WORLD, &status);
    }
 
    /* FOR DEBUGGING PURPOSES ONLY 
     * RUN WITH NPROC = 3 WITH NXY = 16
     * THERE SHOULD BE TWO BLOCKS OF COMPLEMENTARY MATRICES
 
    if (proc_id == 1) {
    printf("\n");
      for (i = 0; i < NXY; i++) {
        for (j = 0; j < cols_per_worker; j++) {
          printf(" %08.0lf", U[0][i][j]);
        }
        printf("\n");
      }
    }
 
    if (proc_id == 2) {
    printf("\n");
      for (i = 0; i < NXY; i++) {
        for (j = cols_per_worker; j < 2*cols_per_worker; j++) {
          printf(" %08.0lf", U[0][i][j]);
        }
        printf("\n");
      }
    }
    */ 
 
 
    /* DETERMINE BORDER ELEMENTS - FIRST & LAST COLUMNS NEED EXTRA CARE:
     * ROW 0 CANNOT EXCHANGE INFORMATION WITH 0-1
     * ROW N CANNOT EXCHANGE INFORMATION WITH N+1
    */
    if (offset == 0) {
      /* DO NOT INCLUDE ZEROTH COLUMN */
      start = 1;
    } else {
      start = offset;
    }
 
    if ((offset + cols_per_worker) == NXY) {
      /* DO NOT INCLUDE THE LAST COLUMN */
      end = offset + cols_per_worker - 2;
    } else {
      end = offset + cols_per_worker - 1;
    }
 
    /* PERFORM TIME STEP [TSTEPS] ITERATIONS.
     * MUST COMMUNICATE BORDER COLUMNS WITH NEIGHBORS.
     * IF IT'S THE FIRST OR THE LAST COLUMN, THEN THE COMMUNICATION TAKES 
     * PLACE ONLY WITH ONE NEIGHBOR
    */
    for(h = 1; h <= TSTEPS; h++) {
      if (lneighbor != NONE) {
        MPI_Send(&U[0][offset], 1, cvector, lneighbor, RNEIGHBOR, MPI_COMM_WORLD);
 
        MPI_Recv(&U[0][offset-1], 1, cvector, lneighbor, LNEIGHBOR, MPI_COMM_WORLD, 
                 &status);
      }
 
      if (rneighbor != NONE) {
        MPI_Send(&U[0][offset+cols_per_worker-1], 1, cvector, rneighbor,
                 LNEIGHBOR, MPI_COMM_WORLD);
 
        MPI_Recv(&U[0][offset+cols_per_worker], 1, cvector, rneighbor,
                 RNEIGHBOR, MPI_COMM_WORLD, &status);
      }
 
      /* UPDATE THE VALUE OF GRID POINTS 
       * REFER TO THE SLIDES FROM WEEK #14 [PART #02]
      */
      for(i = 1; i <= NXY-2; i++) {
        for(j = start; j <= end; j++) {
          Unew[i][j] = U[i][j] + 
                       parameters.cx * (U[i-1][j] + U[i+1][j] - 2 * (U[i][j])) +
                       parameters.cy * (U[i][j-1] + U[i][j+1] - 2 * (U[i][j]));
        }
      }
 
      /* SET U(x,y) TO Unew(x,y) */
      for(i = 1; i <= NXY-2; i++) {
        for(j = start; j <= end; j++) {
          U[i][j] = Unew[i][j];
        }
      }
 
      /* SEND RELEVANT INFORMATION TO MASTER */
      MPI_Send(&offset, 1, MPI_INT, MASTER, DONE, MPI_COMM_WORLD);
      MPI_Send(&cols_per_worker, 1, MPI_INT, MASTER, DONE, MPI_COMM_WORLD);
 
      for(j=0; j < cols_per_worker; j++) {
        MPI_Send(&U[0][offset+j], 1, cvector, MASTER, h*DONE, MPI_COMM_WORLD);
      }
    }
 
  } /* WORKER LOOP ENDS */
 
 
  /* FINALIZE MPI */
  MPI_Finalize();
 
  /* INDICATE THE TERMINATION OF THE PROGRAM */
  return 0;
 
} /* MAIN PROGRAM ENDS */

Program Compilation & Execution

The machine where I am running this calculation, dirac, has 4 processors and has MPICH2 v1.3.1 compiled against GCC v4.1.2 compilers.

[guest@dirac mpi_samples]$ which mpicc
alias mpicc='mpicc -g -Wall -lm'
	~/mpich2/1.3.1/gcc/4.1.2/bin/mpicc
 
[guest@dirac mpi_samples]$ which mpirun
alias mpirun='mpirun -machinefile $HOME/machinefile'
	~/mpich2/1.3.1/gcc/4.1.2/bin/mpirun
 
[guest@dirac mpi_samples]$ mpicc heat_equation_2d.c -o heat_equation_2d.x
 
[guest@dirac mpi_samples]$ mpirun -np 5 ./heat_equation_2d.x
 
  2D Heat Equation [MPI/C]
 
    Number of x-grids        : 20
    Number of y-grids        : 20
    Number of processors     : 5
    Number of workers        : 4
    Number of time steps     : 5000
    Thermal diffusivity (Cx) : 0.050
    Thermal diffusivity (Cy) : 0.050
 
    Temperature distribution (at time t=0)
 
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00032400    00061200    00086400    00108000    00126000    00140400    00151200    00158400    00162000    00162000    00158400    00151200    00140400    00126000    00108000    00086400    00061200    00032400    00000000
    00000000    00061200    00115600    00163200    00204000    00238000    00265200    00285600    00299200    00306000    00306000    00299200    00285600    00265200    00238000    00204000    00163200    00115600    00061200    00000000
    00000000    00086400    00163200    00230400    00288000    00336000    00374400    00403200    00422400    00432000    00432000    00422400    00403200    00374400    00336000    00288000    00230400    00163200    00086400    00000000
    00000000    00108000    00204000    00288000    00360000    00420000    00468000    00504000    00528000    00540000    00540000    00528000    00504000    00468000    00420000    00360000    00288000    00204000    00108000    00000000
    00000000    00126000    00238000    00336000    00420000    00490000    00546000    00588000    00616000    00630000    00630000    00616000    00588000    00546000    00490000    00420000    00336000    00238000    00126000    00000000
    00000000    00140400    00265200    00374400    00468000    00546000    00608400    00655200    00686400    00702000    00702000    00686400    00655200    00608400    00546000    00468000    00374400    00265200    00140400    00000000
    00000000    00151200    00285600    00403200    00504000    00588000    00655200    00705600    00739200    00756000    00756000    00739200    00705600    00655200    00588000    00504000    00403200    00285600    00151200    00000000
    00000000    00158400    00299200    00422400    00528000    00616000    00686400    00739200    00774400    00792000    00792000    00774400    00739200    00686400    00616000    00528000    00422400    00299200    00158400    00000000
    00000000    00162000    00306000    00432000    00540000    00630000    00702000    00756000    00792000    00810000    00810000    00792000    00756000    00702000    00630000    00540000    00432000    00306000    00162000    00000000
    00000000    00162000    00306000    00432000    00540000    00630000    00702000    00756000    00792000    00810000    00810000    00792000    00756000    00702000    00630000    00540000    00432000    00306000    00162000    00000000
    00000000    00158400    00299200    00422400    00528000    00616000    00686400    00739200    00774400    00792000    00792000    00774400    00739200    00686400    00616000    00528000    00422400    00299200    00158400    00000000
    00000000    00151200    00285600    00403200    00504000    00588000    00655200    00705600    00739200    00756000    00756000    00739200    00705600    00655200    00588000    00504000    00403200    00285600    00151200    00000000
    00000000    00140400    00265200    00374400    00468000    00546000    00608400    00655200    00686400    00702000    00702000    00686400    00655200    00608400    00546000    00468000    00374400    00265200    00140400    00000000
    00000000    00126000    00238000    00336000    00420000    00490000    00546000    00588000    00616000    00630000    00630000    00616000    00588000    00546000    00490000    00420000    00336000    00238000    00126000    00000000
    00000000    00108000    00204000    00288000    00360000    00420000    00468000    00504000    00528000    00540000    00540000    00528000    00504000    00468000    00420000    00360000    00288000    00204000    00108000    00000000
    00000000    00086400    00163200    00230400    00288000    00336000    00374400    00403200    00422400    00432000    00432000    00422400    00403200    00374400    00336000    00288000    00230400    00163200    00086400    00000000
    00000000    00061200    00115600    00163200    00204000    00238000    00265200    00285600    00299200    00306000    00306000    00299200    00285600    00265200    00238000    00204000    00163200    00115600    00061200    00000000
    00000000    00032400    00061200    00086400    00108000    00126000    00140400    00151200    00158400    00162000    00162000    00158400    00151200    00140400    00126000    00108000    00086400    00061200    00032400    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
 
 
    Temperature distribution (after 5000 time steps)
 
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000001    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000    00000000
 
 
    Time elapsed (sec)   : 2.5108
 
[guest@dirac mpi_samples]$


Initial and final temperature distribution

Temperature distribution at t=0

Temperature distribution at t=0

Temperature distribution after 5000 time steps

Temperature distribution after 5000 time steps

Temperature distribution

Combined temperature distribution


Animation – using MATLAB

Time lapse of temperature distribution, U(x, y) – 500 time steps, created using the code below.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
% heat_equation_2d.m
% MATLAB 'SCRIPT m' FILE TO READ THE FILE 'heat_equation_2d.dat', 
% ONE DATA SET AT A TIME AND PLOT IT REPEATEDLY TO SIMULATE AN ANIMATION.
%
% TESTED SUCCESSFULLY WITH MATLAB (v7.8.0.347 R2009a; STUDENT VERSION) IN A MACBOOK PRO 
% (OS X 10.6.5) WITH INTEL CORE 2 DUO PROCESSOR (2.4 GHz) & 4GB OF RAM.
%
% FIRST WRITTEN: GOWTHAM; Fri, 10 Dec 2010 23:12:18 -0500
% LAST MODIFIED: GOWTHAM; Sun, 12 Dec 2010 09:34:23 -0500
%
% INPUT:
% 'heat_equation_2d.dat' PRODUCED BY RUNNING 'heat_equation_2d.c'
% AND IT HAS THE FOLLOWING FORMAT:
% 
% TSTEPS=# 
% U(x,y) AT t = 0         [NXY x NXY]
% U(x,y) AT t = 1         [NXY x NXY]
% U(x,y) AT t = 2         [NXY x NXY]
% ...
% U(x,y) AT t = TSTEPS    [NXY x NXY]
%
% FOR AN EXPLANATION OF 'NXY', 'TSTEPS', ETC., PLEASE REFER TO 'heat_equation_2d.c'.
% 
% OUTPUT:
% IF 'mencoder' [AND OTHER REQUIRED UTILITIES] IS AVAILABLE, THIS SCRIP WILL PRODUCE 
% 'heat_equation_2d.avi'. IF NOT, 'New Screen Recording' OPTION MAY BE USED IN
% 'QuickTime Player' TO CREATE THE ANIMATION.
%
 
% CLEAR THE PREVIOUSLY SET VARIABLES
clear
 
% CLEAR THE COMMAND WINDOW
clc
 
% PAUSE FOR 15 SECONDS
% USEFUL WHEN 'mencoder' IS NOT AVAILABLE AS THIS WILL GIVE JUST ENOUGH TIME TO HIDE 
% UNWANTED SCREENS / WINDOWS AND START RECORDING VIA QuickTime Player.
% COMMENT IT OUT IF NOT NEEDED.
pause(15);
 
% OPEN THE DATA FILE FOR READING. 
% SAME BASENAME WILL BE USED FOR THE ANIMATION FILE
filename   = 'heat_equation_2d';
filehandle = fopen(filename '.dat', 'r');
 
% COUNT THE NUMBER OF TIME STEPS - GIVEN IN THE FIRST LINE
TSTEPS = fscanf(filehandle, 'TSTEPS=%d\n', 1);
 
% VARIABLE DECLARATION AND INITIALIZATION
% THIS IS THE NUMBER OF ROWS & COLUMNS IN EACH BLOCK OF DATA
NXY = 20;
 
% x & y GRIDS FOR THE MESH [3D] PLOT
x = 1:1:NXY;
y = 1:1:NXY;
 
% GENERATE X & Y ARRAYS FOR 3D PLOT
% TRANSFORMS THE DOMAIN SPECIFIED BY VECTORS x & y INTO ARRAYS X and Y, 
% WHICH CAN BE USED TO EVALUATE FUNCTIONS OF TWO VARIABLES AND 3D
% MESH/SURFACE PLOTS. THE ROWS OF THE OUTPUT ARRAY X ARE COPIES OF THE
% VECTOR x; COLUMNS OF THE OUTPUT ARRAY Y ARE COPIES OF THE VECTOR y.
% http://www.mathworks.com/help/techdoc/ref/meshgrid.html
[X, Y] = meshgrid(x,y);
 
% LOOP THROUGH THE DATA FILE TO READ 'NXY x NXY' ELEMENTS AT ONCE
% AND IMPORT THEM INTO 'TSTEPS' NUMBER OF MATRICES OF ORDER NXY.
for i = 1:TSTEPS
 
  % INITIALIZE ALL ELEMENTS OF MATRIX [OF ORDER NXY] TO ZERO
  Temperature(i).Z = zeros(NXY, NXY);
 
  % SCAN THE DATA FILE, READ 'NXY x NXY' ELEMENTS, TREAT THEM
  % AS DOUBLE PRECISION AND STORE THEM IN THE MATRIX.
  % SYNTAX:
  % Z = fscanf(filehandle, 'FORMAT', [COLS, ROWS]);
  Temperature(i).Z = fscanf(filehandle, '%f', [NXY, NXY]);
 
  % PLOT THE DATA IN THE ABOVE MATRIX AS A FUNCTION OF X & Y
  figure1 = figure(1);
  mesh(X, Y, Temperature(i).Z)
 
  % HANDLE FOR AXIS PROPERTY SO THAT IT CAN BE MODIFIED TO
  % FIT OUR REQUIREMENS
  ap = gca;
 
  % MODIFY CERTAIN AXIS PROPERTIES APPROPRIATELY
  set(ap, 'XLimMode', 'manual', 'YLimMode', 'manual', ...
          'ZLimMode', 'manual', ...
          'Xlim', [0 20], 'Ylim', [0 20], 'Zlim', [0 1000000], ...
          'FontSize', 14)
 
  % CHANGE THE COLORMAP FROM 'default' TO 'hsv'
  % THE hsv COLORMAP IS CREATED WITH 65536 COLORS
  % http://www.mathworks.com/help/techdoc/ref/colormap.html
  colormap(hsv(65536))
 
  % DISPLAY A COLOR BAR THAT WILL SERVE AS A LEGEND
  colorbar
 
  % AXES LABELS FOR THE PLOT
  xlabel('Grid points along x-axis');
  ylabel('Grid points along y-axis');
  zlabel('U(x, y), Temperature distribution');
 
  % TITLE FOR THE PLOT
  title(['TSTEP= ' num2str(i)]);
 
  % CAPTURE MOVIE FRAME AS AN ARRAY FOR LATER USE.
  % COMMENT THIS LINE IF 'mencoder' IS NOT AVAILABLE.
  % http://www.mathworks.com/help/techdoc/ref/getframe.html
  M(i) = getframe;
 
  % TIME LAG BETWEEN SUCCESSIVE PLOTS
  pause(0.00001);
 
end 
% LOOP THROUGH DATA FILE ENDS
 
% CLOSE THE DATA FILE
fclose(datafilename);
 
 
% CREATE AN ANIMATION USING THE MOVIE FRAMES.
% COMMENT THIS SECTION IF 'mencoder' IS NOT AVAILABLE
 
% CREATE AN AVI [AUDIO/VIDEO INTERLEAVED] FILE FROM THE MOVIE FRAMES ARRAY
% 'None' IS THE ONLY COMPRESSION AVAILABLE OPTION ON UNIX (LIKE) OS.
% http://www.mathworks.com/help/techdoc/ref/movie2avi.html
movie2avi(M, [filename '.avi'], 'compression', 'None', 'fps', 25);
 
% USE 'mencoder' TO COMPRESS THE MOVIE SO THAT IT CAN BE PLAYED ON 
% QuickTime Player OR SOME OTHER SUCH APPLICATION.
% USE eval TO RUN THE MAC/LINUX COMMAND
% http://www.mplayerhq.hu/DOCS/HTML/en/menc-feat-selecting-codec.html
avicompression = ['!mencoder -oac copy -ovc xvid -xvidencopts bitrate=1000 ' filename '.avi -o ' filename '_compressed.avi'];
eval(avicompression);
 
% RENAME THE COMPRESSED AVI
% USE eval TO RUN THE MAC/LINUX COMMAND
avirename = ['!mv ' filename '_compressed.avi ' filename '.avi'];
eval(avirename);