3232// load the PDI header
3333#include <pdi.h>
3434
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
3637int dsize [2 ];
3738
38- /// 2D size of the process grid as [HEIGHT, WIDTH]
39+ // 2D size of the process grid as [HEIGHT, WIDTH]
3940int psize [2 ];
4041
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]
4243int pcoord [2 ];
4344
44- /// the alpha coefficient used in the computation
45+ // the alpha coefficient used in the computation
4546double alpha ;
4647
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
4860 * \param[out] dat the local data to initialize
4961 */
5062void init (double dat [dsize [0 ]][dsize [1 ]])
5163{
5264 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+ }
5486}
5587
5688/** 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]])
6092void iter (double cur [dsize [0 ]][dsize [1 ]], double next [dsize [0 ]][dsize [1 ]])
6193{
6294 int xx , yy ;
63- for (xx = 0 ; xx < dsize [1 ]; ++ xx ) next [0 ][xx ] = cur [0 ][xx ];
6495 for (yy = 1 ; yy < dsize [0 ]- 1 ; ++ yy ) {
65- next [yy ][0 ] = cur [yy ][0 ];
6696 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 ]);
74102 }
75- next [yy ][dsize [1 ]- 1 ] = cur [yy ][dsize [1 ]- 1 ];
76103 }
77- for (xx = 0 ; xx < dsize [1 ]; ++ xx ) next [dsize [0 ]- 1 ][xx ] = cur [dsize [0 ]- 1 ][xx ];
78104}
79105
80- /** Exchanges ghost values with neighbours
106+ /** Exchange ghost values with neighbours
81107 * \param[in] cart_comm the MPI communicator with all processes organized in a 2D Cartesian grid
82108 * \param[in] cur the local data at the current time-step whose ghosts need exchanging
83109 */
@@ -104,8 +130,8 @@ void exchange(MPI_Comm cart_comm, double cur[dsize[0]][dsize[1]])
104130
105131 // send up
106132 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)
109135 cart_comm , & status );
110136
111137 // send to the right
@@ -116,7 +142,7 @@ void exchange(MPI_Comm cart_comm, double cur[dsize[0]][dsize[1]])
116142
117143 // send to the left
118144 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
120146 & cur [1 ][dsize [1 ]- 1 ], 1 , column , rank_source , 100 , // receive last column (ghost)
121147 cart_comm , & status );
122148}
@@ -157,12 +183,12 @@ int main( int argc, char* argv[] )
157183 assert (global_size [1 ]%psize [1 ]== 0 );
158184 assert (psize [1 ]* psize [0 ] == psize_1d );
159185
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)
161187 dsize [0 ] = global_size [0 ]/psize [0 ] + 2 ;
162188 dsize [1 ] = global_size [1 ]/psize [1 ] + 2 ;
163189
164190 // 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 };
166192 MPI_Comm cart_comm ; MPI_Cart_create (main_comm , 2 , psize , cart_period , 1 , & cart_comm );
167193 MPI_Cart_coords (cart_comm , pcoord_1d , 2 , pcoord );
168194
@@ -178,20 +204,20 @@ int main( int argc, char* argv[] )
178204 int ii = 0 ;
179205
180206 // 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+ //...
183209 PDI_share ("pcoord" , pcoord , PDI_OUT );
184210 PDI_reclaim ("pcoord" );
185211 PDI_share ("dsize" , dsize , PDI_OUT );
186212 PDI_reclaim ("dsize" );
187213 PDI_share ("psize" , psize , PDI_OUT );
188214 PDI_reclaim ("psize" );
189- PDI_share ("main_field" , cur , PDI_OUT );
190- PDI_reclaim ("main_field" );
191215
192216 // the main loop
193- for (; ii < 10 ; ++ ii ) {
217+ for (; ii < 4 ; ++ ii ) {
194218 // share the loop counter & main field at each iteration
219+ //*** use PDI_multi_expose to replace PDI_share + event + reclaim
220+ //...
195221 PDI_share ("ii" , & ii , PDI_OUT );
196222 PDI_share ("main_field" , cur , PDI_OUT );
197223 PDI_event ("loop" );
@@ -208,11 +234,13 @@ int main( int argc, char* argv[] )
208234 double (* tmp )[dsize [1 ]] = cur ; cur = next ; next = tmp ;
209235 }
210236 // 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+ //...
212239 PDI_share ("ii" , & ii , PDI_OUT );
240+ PDI_share ("main_field" , cur , PDI_OUT );
213241 PDI_event ("finalization" );
214- PDI_reclaim ("ii" );
215242 PDI_reclaim ("main_field" );
243+ PDI_reclaim ("ii" );
216244
217245 // finalize PDI
218246 PDI_finalize ();
@@ -230,3 +258,4 @@ int main( int argc, char* argv[] )
230258 fprintf (stderr , "[%d] SUCCESS\n" , pcoord_1d );
231259 return EXIT_SUCCESS ;
232260}
261+
0 commit comments