3333// load the PDI header
3434#include <pdi.h>
3535
36- /// size of the local data as [HEIGHT, WIDTH] including ghosts & boundary constants
36+ // size of the local data as [HEIGHT, WIDTH] including the number of ghost layers
37+ // for communications or boundary conditions
3738int dsize [2 ];
3839
39- /// 2D size of the process grid as [HEIGHT, WIDTH]
40+ // 2D size of the process grid as [HEIGHT, WIDTH]
4041int psize [2 ];
4142
42- /// 2D rank of the local process in the process grid as [YY, XX]
43+ // 2D rank of the local process in the process grid as [YY, XX]
4344int pcoord [2 ];
4445
45- /// the alpha coefficient used in the computation
46+ // the alpha coefficient used in the computation
4647double alpha ;
4748
4849double L = 1.0 ;
50+ // definition of the source
51+ // the source corresponds to a disk of an uniform value
52+ // source1: center=(0.4,0.4), radius=0.2 and value=100
4953double source1 [4 ]= {0.4 , 0.4 , 0.2 , 100 };
54+ // source2: center=(0.8,0.7), radius=0.1 and value=200
5055double source2 [4 ]= {0.7 , 0.8 , 0.1 , 200 };
56+ // the order of the coordinates of the center (XX,YY) is inverted in the vector
5157
5258FILE * pFile2 = NULL ;
5359
@@ -77,7 +83,9 @@ void close_file(void)
7783 fclose (pFile2 );
7884}
7985
80- /** Initialize the data all to 0 except for the left border (XX==0) initialized to 1 million
86+ /** Initialize all the data to 0, with the exception of each cells
87+ * whose center (cpos_x,cpos_y) is inside of the disks
88+ * defined by source1 or source2
8189 * \param[out] dat the local data to initialize
8290 */
8391void init (double dat [dsize [0 ]][dsize [1 ]])
@@ -87,14 +95,19 @@ void init(double dat[dsize[0]][dsize[1]])
8795 double dx = L / ((dsize [1 ]- 2 ) * psize [1 ]) ;
8896
8997 double cpos_x ,cpos_y ;
98+ double square_dist1 , square_dist2 ;
9099 for (int yy = 0 ; yy < dsize [0 ];++ yy ) {
91100 cpos_y = (yy + pcoord [0 ]* (dsize [0 ]- 2 ))* dy - 0.5 * dy ;
92101 for (int xx = 0 ; xx < dsize [1 ];++ xx ) {
93102 cpos_x = (xx + pcoord [1 ]* (dsize [1 ]- 2 ))* dx - 0.5 * dx ;
94- if ((cpos_y - source1 [0 ])* (cpos_y - source1 [0 ]) + (cpos_x - source1 [1 ])* (cpos_x - source1 [1 ]) <= source1 [2 ]* source1 [2 ]) {
103+ square_dist1 = ( cpos_y - source1 [0 ] ) * ( cpos_y - source1 [0 ] )
104+ + ( cpos_x - source1 [1 ] ) * ( cpos_x - source1 [1 ] );
105+ if (square_dist1 <= source1 [2 ] * source1 [2 ]) {
95106 dat [yy ][xx ] = source1 [3 ];
96107 }
97- if ((cpos_y - source2 [0 ])* (cpos_y - source2 [0 ]) + (cpos_x - source2 [1 ])* (cpos_x - source2 [1 ]) <= source2 [2 ]* source2 [2 ]) {
108+ square_dist2 = ( cpos_y - source2 [0 ] ) * ( cpos_y - source2 [0 ] )
109+ + ( cpos_x - source2 [1 ] ) * ( cpos_x - source2 [1 ] );
110+ if (square_dist2 <= source2 [2 ] * source2 [2 ]) {
98111 dat [yy ][xx ] = source2 [3 ];
99112 }
100113 }
@@ -119,7 +132,7 @@ void iter(double cur[dsize[0]][dsize[1]], double next[dsize[0]][dsize[1]])
119132 }
120133}
121134
122- /** Exchanges ghost values with neighbours
135+ /** Exchange ghost values with neighbours
123136 * \param[in] cart_comm the MPI communicator with all processes organized in a 2D Cartesian grid
124137 * \param[in] cur the local data at the current time-step whose ghosts need exchanging
125138 */
@@ -146,8 +159,8 @@ void exchange(MPI_Comm cart_comm, double cur[dsize[0]][dsize[1]])
146159
147160 // send up
148161 MPI_Cart_shift (cart_comm , 0 , -1 , & rank_source , & rank_dest );
149- MPI_Sendrecv (& cur [1 ][1 ], 1 , row , rank_dest , 100 , // send column after ghost
150- & cur [dsize [0 ]- 1 ][1 ], 1 , row , rank_source , 100 , // receive last column (ghost)
162+ MPI_Sendrecv (& cur [1 ][1 ], 1 , row , rank_dest , 100 , // send row after ghost
163+ & cur [dsize [0 ]- 1 ][1 ], 1 , row , rank_source , 100 , // receive last row (ghost)
151164 cart_comm , & status );
152165
153166 // send to the right
@@ -158,17 +171,14 @@ void exchange(MPI_Comm cart_comm, double cur[dsize[0]][dsize[1]])
158171
159172 // send to the left
160173 MPI_Cart_shift (cart_comm , 1 , -1 , & rank_source , & rank_dest );
161- MPI_Sendrecv (& cur [1 ][1 ], 1 , column , rank_dest , 100 , // send column after ghost
174+ MPI_Sendrecv (& cur [1 ][1 ], 1 , column , rank_dest , 100 , // send column after ghost
162175 & cur [1 ][dsize [1 ]- 1 ], 1 , column , rank_source , 100 , // receive last column (ghost)
163176 cart_comm , & status );
164177}
165178
166179int main ( int argc , char * argv [] )
167180{
168181 MPI_Init (& argc , & argv );
169- srand ( time ( NULL ) );
170-
171- int switch_iter_value [10 ] = {20 , 35 , 50 , 55 , 60 , 35 , 25 , 20 , 15 , 60 };
172182
173183 // load the configuration tree
174184 PC_tree_t conf = PC_parse_path ("ex12.yml" );
@@ -202,7 +212,7 @@ int main( int argc, char* argv[] )
202212 assert (global_size [1 ]%psize [1 ]== 0 );
203213 assert (psize [1 ]* psize [0 ] == psize_1d );
204214
205- // compute the local data-size with space for ghosts and boundary constants
215+ // compute the local data-size (the number of ghost layers is 2 for each coordinate)
206216 dsize [0 ] = global_size [0 ]/psize [0 ] + 2 ;
207217 dsize [1 ] = global_size [1 ]/psize [1 ] + 2 ;
208218
0 commit comments