32
32
// load the PDI header
33
33
#include <pdi.h>
34
34
35
- /// size of the local data as [HEIGHT, WIDTH] including ghosts & boundary constants
35
+ // size of the local data as [HEIGHT, WIDTH] including the number of ghost layers
36
+ // for communications or boundary conditions
36
37
int dsize [2 ];
37
38
38
- /// 2D size of the process grid as [HEIGHT, WIDTH]
39
+ // 2D size of the process grid as [HEIGHT, WIDTH]
39
40
int psize [2 ];
40
41
41
- /// 2D rank of the local process in the process grid as [YY, XX]
42
+ // 2D rank of the local process in the process grid as [YY, XX]
42
43
int pcoord [2 ];
43
44
44
- /// the alpha coefficient used in the computation
45
+ // the alpha coefficient used in the computation
45
46
double alpha ;
46
47
47
- /** Initialize the data all to 0 except for the left border (XX==0) initialized to 1 million
48
+ double L = 1.0 ;
49
+ // definition of the source
50
+ // the source corresponds to a disk of an uniform value
51
+ // source1: center=(0.4,0.4), radius=0.2 and value=100
52
+ double source1 [4 ]= {0.4 , 0.4 , 0.2 , 100 };
53
+ // source2: center=(0.8,0.7), radius=0.1 and value=200
54
+ double source2 [4 ]= {0.7 , 0.8 , 0.1 , 200 };
55
+ // the order of the coordinates of the center (XX,YY) is inverted in the vector
56
+
57
+ /** Initialize all the data to 0, with the exception of a given cell
58
+ * whose center (cpos_x,cpos_y) is inside of the disks
59
+ * defined by source1 or source2
48
60
* \param[out] dat the local data to initialize
49
61
*/
50
62
void init (double dat [dsize [0 ]][dsize [1 ]])
51
63
{
52
64
for (int yy = 0 ; yy < dsize [0 ]; ++ yy ) for (int xx = 0 ; xx < dsize [1 ]; ++ xx ) dat [yy ][xx ] = 0 ;
53
- if ( pcoord [1 ] == 0 ) for (int yy = 0 ; yy < dsize [0 ]; ++ yy ) dat [yy ][0 ] = 1000000 ;
65
+ double dy = L / ((dsize [0 ]- 2 ) * psize [0 ]) ;
66
+ double dx = L / ((dsize [1 ]- 2 ) * psize [1 ]) ;
67
+
68
+ double cpos_x ,cpos_y ;
69
+ double square_dist1 , square_dist2 ;
70
+ for (int yy = 0 ; yy < dsize [0 ];++ yy ) {
71
+ cpos_y = (yy + pcoord [0 ]* (dsize [0 ]- 2 ))* dy - 0.5 * dy ;
72
+ for (int xx = 0 ; xx < dsize [1 ];++ xx ) {
73
+ cpos_x = (xx + pcoord [1 ]* (dsize [1 ]- 2 ))* dx - 0.5 * dx ;
74
+ square_dist1 = ( cpos_y - source1 [0 ] ) * ( cpos_y - source1 [0 ] )
75
+ + ( cpos_x - source1 [1 ] ) * ( cpos_x - source1 [1 ] );
76
+ if (square_dist1 <= source1 [2 ] * source1 [2 ]) {
77
+ dat [yy ][xx ] = source1 [3 ];
78
+ }
79
+ square_dist2 = ( cpos_y - source2 [0 ] ) * ( cpos_y - source2 [0 ] )
80
+ + ( cpos_x - source2 [1 ] ) * ( cpos_x - source2 [1 ] );
81
+ if (square_dist2 <= source2 [2 ] * source2 [2 ]) {
82
+ dat [yy ][xx ] = source2 [3 ];
83
+ }
84
+ }
85
+ }
54
86
}
55
87
56
88
/** Compute the values at the next time-step based on the values at the current time-step
@@ -60,24 +92,18 @@ void init(double dat[dsize[0]][dsize[1]])
60
92
void iter (double cur [dsize [0 ]][dsize [1 ]], double next [dsize [0 ]][dsize [1 ]])
61
93
{
62
94
int xx , yy ;
63
- for (xx = 0 ; xx < dsize [1 ]; ++ xx ) next [0 ][xx ] = cur [0 ][xx ];
64
95
for (yy = 1 ; yy < dsize [0 ]- 1 ; ++ yy ) {
65
- next [yy ][0 ] = cur [yy ][0 ];
66
96
for (xx = 1 ; xx < dsize [1 ]- 1 ; ++ xx ) {
67
- next [yy ][xx ] =
68
- (1. - 4. * alpha ) * cur [yy ][xx ]
69
- + alpha * ( cur [yy ][xx - 1 ]
70
- + cur [yy ][xx + 1 ]
71
- + cur [yy - 1 ][xx ]
72
- + cur [yy + 1 ][xx ]
73
- );
97
+ next [yy ][xx ] = (1. - 4. * alpha ) * cur [yy ][xx ]
98
+ + alpha * ( cur [yy ][xx - 1 ]
99
+ + cur [yy ][xx + 1 ]
100
+ + cur [yy - 1 ][xx ]
101
+ + cur [yy + 1 ][xx ]);
74
102
}
75
- next [yy ][dsize [1 ]- 1 ] = cur [yy ][dsize [1 ]- 1 ];
76
103
}
77
- for (xx = 0 ; xx < dsize [1 ]; ++ xx ) next [dsize [0 ]- 1 ][xx ] = cur [dsize [0 ]- 1 ][xx ];
78
104
}
79
105
80
- /** Exchanges ghost values with neighbours
106
+ /** Exchange ghost values with neighbours
81
107
* \param[in] cart_comm the MPI communicator with all processes organized in a 2D Cartesian grid
82
108
* \param[in] cur the local data at the current time-step whose ghosts need exchanging
83
109
*/
@@ -104,8 +130,8 @@ void exchange(MPI_Comm cart_comm, double cur[dsize[0]][dsize[1]])
104
130
105
131
// send up
106
132
MPI_Cart_shift (cart_comm , 0 , -1 , & rank_source , & rank_dest );
107
- MPI_Sendrecv (& cur [1 ][1 ], 1 , row , rank_dest , 100 , // send column after ghost
108
- & cur [dsize [0 ]- 1 ][1 ], 1 , row , rank_source , 100 , // receive last column (ghost)
133
+ MPI_Sendrecv (& cur [1 ][1 ], 1 , row , rank_dest , 100 , // send row after ghost
134
+ & cur [dsize [0 ]- 1 ][1 ], 1 , row , rank_source , 100 , // receive last row (ghost)
109
135
cart_comm , & status );
110
136
111
137
// send to the right
@@ -116,7 +142,7 @@ void exchange(MPI_Comm cart_comm, double cur[dsize[0]][dsize[1]])
116
142
117
143
// send to the left
118
144
MPI_Cart_shift (cart_comm , 1 , -1 , & rank_source , & rank_dest );
119
- MPI_Sendrecv (& cur [1 ][1 ], 1 , column , rank_dest , 100 , // send column after ghost
145
+ MPI_Sendrecv (& cur [1 ][1 ], 1 , column , rank_dest , 100 , // send column after ghost
120
146
& cur [1 ][dsize [1 ]- 1 ], 1 , column , rank_source , 100 , // receive last column (ghost)
121
147
cart_comm , & status );
122
148
}
@@ -157,12 +183,12 @@ int main( int argc, char* argv[] )
157
183
assert (global_size [1 ]%psize [1 ]== 0 );
158
184
assert (psize [1 ]* psize [0 ] == psize_1d );
159
185
160
- // compute the local data-size with space for ghosts and boundary constants
186
+ // compute the local data-size (the number of ghost layers is 2 for each coordinate)
161
187
dsize [0 ] = global_size [0 ]/psize [0 ] + 2 ;
162
188
dsize [1 ] = global_size [1 ]/psize [1 ] + 2 ;
163
189
164
190
// create a 2D Cartesian MPI communicator & get our coordinate (rank) in it
165
- int cart_period [2 ] = { 0 , 0 };
191
+ int cart_period [2 ] = { 1 , 1 };
166
192
MPI_Comm cart_comm ; MPI_Cart_create (main_comm , 2 , psize , cart_period , 1 , & cart_comm );
167
193
MPI_Cart_coords (cart_comm , pcoord_1d , 2 , pcoord );
168
194
@@ -178,20 +204,20 @@ int main( int argc, char* argv[] )
178
204
int ii = 0 ;
179
205
180
206
// share useful configuration bits with PDI
181
- PDI_share ( "ii" , & ii , PDI_OUT );
182
- PDI_reclaim ( "ii" );
207
+ //*** use PDI_expose to replace PDI_share + reclaim
208
+ //...
183
209
PDI_share ("pcoord" , pcoord , PDI_OUT );
184
210
PDI_reclaim ("pcoord" );
185
211
PDI_share ("dsize" , dsize , PDI_OUT );
186
212
PDI_reclaim ("dsize" );
187
213
PDI_share ("psize" , psize , PDI_OUT );
188
214
PDI_reclaim ("psize" );
189
- PDI_share ("main_field" , cur , PDI_OUT );
190
- PDI_reclaim ("main_field" );
191
215
192
216
// the main loop
193
- for (; ii < 10 ; ++ ii ) {
217
+ for (; ii < 4 ; ++ ii ) {
194
218
// share the loop counter & main field at each iteration
219
+ //*** use PDI_multi_expose to replace PDI_share + event + reclaim
220
+ //...
195
221
PDI_share ("ii" , & ii , PDI_OUT );
196
222
PDI_share ("main_field" , cur , PDI_OUT );
197
223
PDI_event ("loop" );
@@ -208,11 +234,13 @@ int main( int argc, char* argv[] )
208
234
double (* tmp )[dsize [1 ]] = cur ; cur = next ; next = tmp ;
209
235
}
210
236
// finally share the main field as well as the loop counter after the loop
211
- PDI_share ("main_field" , cur , PDI_OUT );
237
+ //*** use PDI_multi_expose to replace PDI_share + event + reclaim
238
+ //...
212
239
PDI_share ("ii" , & ii , PDI_OUT );
240
+ PDI_share ("main_field" , cur , PDI_OUT );
213
241
PDI_event ("finalization" );
214
- PDI_reclaim ("ii" );
215
242
PDI_reclaim ("main_field" );
243
+ PDI_reclaim ("ii" );
216
244
217
245
// finalize PDI
218
246
PDI_finalize ();
@@ -230,3 +258,4 @@ int main( int argc, char* argv[] )
230
258
fprintf (stderr , "[%d] SUCCESS\n" , pcoord_1d );
231
259
return EXIT_SUCCESS ;
232
260
}
261
+
0 commit comments