Newer
Older
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
/* Setup routines for heat equation solver */
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <assert.h>
#include <mpi.h>
#include "heat.h"
#include "pngwriter.h"
#define NSTEPS 500 // Default number of iteration steps
/* Initialize the heat equation solver */
void initialize(int argc, char *argv[], field *current,
field *previous, int *nsteps, parallel_data *parallel,
int *iter0)
{
/*
* Following combinations of command line arguments are possible:
* No arguments: use default field dimensions and number of time steps
* One argument: read initial field from a given file
* Two arguments: initial field from file and number of time steps
* Three arguments: field dimensions (rows,cols) and number of time steps
*/
int rows = 2000; //!< Field dimensions with default values
int cols = 2000;
char input_file[64]; //!< Name of the optional input file
int read_file = 0;
*nsteps = NSTEPS;
*iter0 = 0;
switch (argc) {
case 1:
/* Use default values */
break;
case 2:
/* Read initial field from a file */
strncpy(input_file, argv[1], 64);
read_file = 1;
break;
case 3:
/* Read initial field from a file */
strncpy(input_file, argv[1], 64);
read_file = 1;
/* Number of time steps */
*nsteps = atoi(argv[2]);
break;
case 4:
/* Field dimensions */
rows = atoi(argv[1]);
cols = atoi(argv[2]);
/* Number of time steps */
*nsteps = atoi(argv[3]);
break;
default:
printf("Unsupported number of command line arguments\n");
exit(-1);
}
// Check if checkpoint exists
if (!access(CHECKPOINT, F_OK)) {
read_restart(current, parallel, iter0);
set_field_dimensions(previous, current->nx_full, current->ny_full,
parallel);
allocate_field(previous);
if (parallel->rank == 0)
printf("Restarting from an earlier checkpoint saved"
" at iteration %d.\n", *iter0);
copy_field(current, previous);
} else if (read_file) {
read_field(current, previous, input_file, parallel);
} else {
parallel_setup(parallel, rows, cols);
set_field_dimensions(current, rows, cols, parallel);
set_field_dimensions(previous, rows, cols, parallel);
generate_field(current, parallel);
allocate_field(previous);
copy_field(current, previous);
}
}
/* Generate initial temperature field. Pattern is disc with a radius
* of nx_full / 6 in the center of the grid.
* Boundary conditions are (different) constant temperatures outside the grid */
void generate_field(field *temperature, parallel_data *parallel)
{
int i, j, ind, width;
double radius;
int dx, dy;
int dims[2], coords[2], periods[2];
/* Allocate the temperature array, note that
* we have to allocate also the ghost layers */
temperature->data =
malloc_2d(temperature->nx + 2, temperature->ny + 2);
MPI_Cart_get(parallel->comm, 2, dims, periods, coords);
/* Radius of the source disc */
radius = temperature->nx_full / 6.0;
width = temperature->ny + 2;
for (i = 0; i < temperature->nx + 2; i++) {
for (j = 0; j < temperature->ny + 2; j++) {
/* Distance of point i, j from the origin */
dx = i + coords[0] * temperature->nx -
temperature->nx_full / 2 + 1;
dy = j + coords[1] * temperature->ny -
temperature->ny_full / 2 + 1;
ind = idx(i, j, width);
if (dx * dx + dy * dy < radius * radius) {
temperature->data[ind] = 5.0;
} else {
temperature->data[ind] = 65.0;
}
}
}
/* Boundary conditions */
// Left boundary
if (coords[1] == 0) {
for (i = 0; i < temperature->nx + 2; i++) {
ind = idx(i, 0, width);
temperature->data[ind] = 20.0;
}
}
// Right boundary
if (coords[1] == dims[1] - 1) {
for (i = 0; i < temperature->nx + 2; i++) {
ind = idx(i, temperature->ny + 1, width);
temperature->data[ind] = 70.0;
}
}
// Upper boundary
if (coords[0] == 0) {
for (j = 0; j < temperature->ny + 2; j++) {
ind = idx(0, j, width);
temperature->data[ind] = 85.0;
}
}
// Lower boundary
if (coords[0] == dims[0] - 1) {
for (j = 0; j < temperature->ny + 2; j++) {
ind = idx(temperature->nx + 1, j, width);
temperature->data[ind] = 5.0;
}
}
}
/* Set dimensions of the field. Note that the nx is the size of the first
* dimension and ny the second. */
void set_field_dimensions(field *temperature, int nx, int ny,
parallel_data *parallel)
{
int nx_local, ny_local;
int dims[2], coords[2], periods[2];
MPI_Cart_get(parallel->comm, 2, dims, periods, coords);
nx_local = nx / dims[0];
ny_local = ny / dims[1];
temperature->dx = DX;
temperature->dy = DY;
temperature->nx = nx_local;
temperature->ny = ny_local;
temperature->nx_full = nx;
temperature->ny_full = ny;
}
void parallel_setup(parallel_data *parallel, int nx, int ny)
{
int nx_local;
int ny_local;
int world_size;
int dims[2] = {0, 0};
int periods[2] = { 0, 0 };
/* Set grid dimensions */
MPI_Comm_size(MPI_COMM_WORLD, &world_size);
MPI_Dims_create(world_size, 2, dims);
nx_local = nx / dims[0];
ny_local = ny / dims[1];
if (nx_local * dims[0] != nx) {
printf("Cannot divide grid evenly to processors in x-direction "
"%d x %d != %d\n", nx_local, dims[0], nx);
MPI_Abort(MPI_COMM_WORLD, -2);
}
if (ny_local * dims[1] != ny) {
printf("Cannot divide grid evenly to processors in y-direction "
"%d x %d != %d\n", ny_local, dims[1], ny);
MPI_Abort(MPI_COMM_WORLD, -2);
}
/* Create cartesian communicator */
MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, 1, ¶llel->comm);
MPI_Cart_shift(parallel->comm, 0, 1, ¶llel->nup, ¶llel->ndown);
MPI_Cart_shift(parallel->comm, 1, 1, ¶llel->nleft,
¶llel->nright);
MPI_Comm_size(parallel->comm, ¶llel->size);
MPI_Comm_rank(parallel->comm, ¶llel->rank);
if (parallel->rank == 0) {
printf("Using domain decomposition %d x %d\n", dims[0], dims[1]);
printf("Local domain size %d x %d\n", nx_local, ny_local);
}
/* Create datatypes for halo exchange */
MPI_Type_vector(nx_local + 2, 1, ny_local + 2, MPI_DOUBLE,
¶llel->columntype);
MPI_Type_contiguous(ny_local + 2, MPI_DOUBLE, ¶llel->rowtype);
MPI_Type_commit(¶llel->columntype);
MPI_Type_commit(¶llel->rowtype);
/* Create datatype for subblock needed in text I/O
* Rank 0 uses datatype for receiving data into full array while
* other ranks use datatype for sending the inner part of array */
int sizes[2] = {nx_local + 2, ny_local + 2};
int subsizes[2] = { nx_local, ny_local };
int offsets[2] = {1, 1};
if (parallel->rank == 0) {
sizes[0] = nx;
sizes[1] = ny;
offsets[0] = 0;
offsets[1] = 0;
}
MPI_Type_create_subarray(2, sizes, subsizes, offsets, MPI_ORDER_C,
MPI_DOUBLE, ¶llel->subarraytype);
MPI_Type_commit(¶llel->subarraytype);
/* Create datatypes for restart I/O
* For boundary ranks also the ghost layer (boundary condition)
* is written */
int coords[2];
MPI_Cart_coords(parallel->comm, parallel->rank, 2, coords);
sizes[0] = nx + 2;
sizes[1] = ny + 2;
offsets[0] = 1 + coords[0] * nx_local;
offsets[1] = 1 + coords[1] * ny_local;
if (coords[0] == 0) {
offsets[0] -= 1;
subsizes[0] += 1;
}
if (coords[0] == dims[0] - 1) {
subsizes[0] += 1;
}
if (coords[1] == 0) {
offsets[1] -= 1;
subsizes[1] += 1;
}
if (coords[1] == dims[1] - 1) {
subsizes[1] += 1;
}
MPI_Type_create_subarray(2, sizes, subsizes, offsets, MPI_ORDER_C,
MPI_DOUBLE, ¶llel->filetype);
MPI_Type_commit(¶llel->filetype);
sizes[0] = nx_local + 2;
sizes[1] = ny_local + 2;
offsets[0] = 1;
offsets[1] = 1;
if (coords[0] == 0) {
}
if (coords[1] == 0) {
}
MPI_Type_create_subarray(2, sizes, subsizes, offsets, MPI_ORDER_C,
MPI_DOUBLE, ¶llel->restarttype);
MPI_Type_commit(¶llel->restarttype);
}
/* Deallocate the 2D arrays of temperature fields */
void finalize(field *temperature1, field *temperature2,
parallel_data *parallel)
{
free_2d(temperature1->data);
free_2d(temperature2->data);
MPI_Type_free(¶llel->rowtype);
MPI_Type_free(¶llel->columntype);
MPI_Type_free(¶llel->subarraytype);
MPI_Type_free(¶llel->restarttype);
MPI_Type_free(¶llel->filetype);
}